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Abstract 


This  report  summarizes  recent  work  directed  toward  the  development  of  an  effective 
finite  element-based  eigensolution  capability  for  coupled  fluid /structure  interaction  prob¬ 
lems  in  which  the  fluid  compressibility  is  included.  Solution  of  these  coupled  problems  is 
essential  in  the  evaluation  of  structural  acoustic  transmission  characteristics  and  enclosure 
acoustics.  The  methodologies  are  also  of  direct  relevance  to  the  acoustic  radiation  problem 
as  well  as  other  forms  of  time  domain  analysis.  A  comprehensive  review  of  the  four  general 
methods  of  treating  this  problem  is  presented  with  the  objective  of  establishing  future  di¬ 
rections  most  compatible  with  existing  finite  element  analysis  capability.  Several  example 
applications  of  the  different  methods  are  presented  to  highlight  practical  considerations. 
Of  the  formulations  reviewed,  a  three-field  system  appears  to  have  the  most  potential  as 
a  general  purpose  approach,  although  considerable  effort  is  required  in  the  development  of 
a  suitable  eigensolver  for  the  system  of  choice.  The  report  also  presents  a  discussion  of 
a  number  of  methods  of  incorporating  concentrated  damping  in  finite  element  models  of 
acoustic  domains. 


Resume 


Le  present  rapport  resume  les  travaux  effectuds  pour  mettre  au  point  une  solution 
caractdristique  aux  elements  finis  efficace  aux  problfcmes  d'interaction  fluide/structure 
couples,  solution  qui  tient  compte  de  la  compressibility  du  fluide.  D  est  essentiel  de 
rdsoudre  ces  probtemes  de  systfemes  couples  pour  ^valuer  les  caractdristiques  de 
transmission  de  sons  par  les  structures  et  1'acoustique  d'une  enceinte.  Les  m6thodes  sont 
en  outre  directement  lides  au  probleme  de  rayonnement  acoustique  ainsi  qu'ii  d'autres 
formes  d'analyse  du  domaine  temporel.  Une  etude  compr6hensif  des  quatre  mdthodes 
generates  de  traitement  de  ce  probleme  est  presentee  dans  le  but  d'dtablirdes  orientations 
futures  qui  soient  des  plus  compatibles  avec  les  moyens  actuels  d'analyse  par  elements 
finis.  Plusieurs  applications  des  differentes  methodes  sont  donn6es  en  exemple  pour  faire 
ressortir  des  facteurs  d'ordre  pratique.  Des  formulations  presentees,  un  systfcme  k  trois 
champs  semble  le  plus  prometteur  comme  mdthode  de  resolution  generate,  meme  si  un 
effort  considerable  est  necessaire  pour  trouver  une  solution  caracteristique  appropriee  au 
systeme  choisi.  Le  rapport  renferme  aussi  une  analyse  de  plusieurs  methodes  pour  integrer 
des  amortissemems  concentres  dans  des  modeies  de  domaines  acousdques  aux  elements 
finis. 
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Notation 

In  the  presentation  of  many  different  formulations,  consistency  of  notation  is  difficult. 
In  this  report,  every  endeavour  has  been  made  to  retain  the  same  symbol  for  the  same  type 
of  discrete  operation,  although  the  terms  involved  may  be  quite  different.  Thus  for  example, 
L  always  refers  to  a  coupling  matrix  generated  from  cross-products  of  structural  and  fluid 
variable  representations,  whether  based  on  modal  or  interpolation  functions.  Vectors  are 
printed  in  boldface  type. 

A  Admittance  function,  boundary  area 

Ay ,  A-i ,  As  Coupling  matrices 
B,  By  Fluid  coupling  matrices 

c,  ca  Sound  speed  in  fluid  medium,  absorptive  material 

C  Absorptive  boundary  damping  matrix 

D  Elasticity  matrix 

D,  Distributed  structural  damping  matrix 

Df  Distributed  fluid  damping  matrix 

E  Young's  modulus 

F  Load  vector 

g  gravitational  acceleration 

H,  Hy  Fluid  stiffness  type  matrices 

Hr  Rotational  stiffness  matrix 

I  Identity  matrix 

j  y/=I 

K  Structural  stiffness  matrix 

K,  Structural  factor  of  absorptive  material 

K  Generalized  stiffness  <pTK<p 

L,  Ly,Li  Fluid/structure  coupling  matrix 


IV 


M 

M 

Ma 

n 

N 

P,Po 

Pn 

Q 

Qa 
R 
u 
t n 

V 

X 

Z 

a 

P 

i  r#,  r  v 

6 

? 

V 
A 
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Structural  mass  matrix 

Generalized  mass  matrix  <pTM<p 

Added  mass  matrix 

Unit  normal  vector 

Finite  element  interpolation  function 

Dynamic,  hydrostatic  pressure 

Non-dimensional  pressure  coefficient 

Fluid  mass  type  matrix 

Mass  type  matrix  for  absorptive  material 

Absorptive  material  flow  resistivity 

Structural  displacement  vector 

Displacement  component 

Fluid  particle,  boundary  velocity 

Body  force  vector 

Impedance  function 

Penalty  function  coefficient 

Fluid  bulk  modulus 

Wetted,  absorptive,  specified  velocity  boundary 
Variational  operator 

Generalized  structural  displacement  coordinates 
Generalized  fluid  displacement  coordinates 
Laml  constant 

Lam£  constant  (shear  modulus) 

Poisson’s  ratio 
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p,p$ 

4> 

<P;<Pf 

u 

ua,ut 

n 

Clf  ,Qt 
* 

V.,V 


Dynamic,  static,  fluid  density 
Velocity  potential  function 
Structural,  fluid  displacement  mode 
Displacement  potential  function 
Coupled  natural  frequency 

Rigid-walled  fluid  cavity,  structural  natural  frequency 

Absorptive  material  porosity 

Fluid,  structural  domain 

Cross-product  operator 

Divergence,  gradient  operator 
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1  Introduction 


This  report  discusses  finite  element  methodologies  for  the  prediction  of  natural  modes  of 
vibration  of  coupled  fluid/structure  systems  in  which  the  fluid  compressibility  is  included. 
These  formulations  find  applicability  in  a  wide  variety  of  important  problems,  such  as  the 
evaluation  of  acoustic  transmission  characteristics,  prediction  of  enclosure  acoustics,  and 
vibration  prediction  for  submerged  structures.  The  solution  of  the  coupled  eigenproblem  is 
also  often  of  direct  relevance  to  other  forms  of  steady  state  and  transient  response  analysis. 

The  relatively  straightforward  application  of  finite  element  methods  to  the  acoustic 
modal  analysis  of  bounded  domains  and  the  extent  to  which  this  method  has  become  the 
standard  for  structural  analysis  raise  few  suspicions  that  the  treatment  of  coupled  compress¬ 
ible  fluid/structure  systems  would  be  problematic.  That  this  is  not  the  case  is  apparent 
from  the  continued  discussion  of  this  topic  in  the  literature,  a  full  twenty  years  after  the  ini¬ 
tial  approaches  were  undertaken.  While  the  original  intent  of  this  investigation  was  simply 
to  establish  a  viable  approach  to  the  coupled  eigenvalue  problem  within  the  framework  of 
an  existing  structural  finite  element  analysis  program  suite,  the  large  number  of  seemingly 
disparate  approaches  proposed  in  the  literature  has  led  to  a  much  more  in-depth  review  of 
this  problem. 

Finite  element  formulation  of  an  initially  undamped  coupled  eigenvalue  problem  leads  to 
one  or  more  of  the  following  problems:  the  presence  of  unsymmetric  coefficient  matrices;  the 
presence  of  a  pseudo  damping  matrix;  inefficiency  arising  from  the  inclusion  of  extra  vari¬ 
ables;  highly  singular  and/ or  indefinite  coefficient  matrices  or  the  presence  of  numerous  zero 
strain  energy  modes;  necessity  of  matrix  inversion;  large  bandwidths  and  ill-conditioning. 
The  occurrence  of  most  of  these  characteristics  is  a  function  of  the  particular  choice  of  the 
variable  or  variables  used  to  describe  the  behaviour  in  the  individual  structural  and  fluid 
domains.  To  exploit  conventional  structural  finite  element  solution  systems  to  advantage, 
displacements  are  invariably  chosen  for  the  structural  domain,  although  it  could  be  argued 
that  force  variables  would  be  a  more  logical  choice.  In  the  fluid  region,  one  finds  formu¬ 
lations  based  on  one  or  more  of  the  fundamental  variables  of  displacement,  displacement 
potential,  pressure,  and  velocity  potential.  It  is  of  note  that  formulations  which  appear  to 
be  quite  different  can  often  be  derived  from  each  other  using  the  simple  relationships  which 
must  exist  between  these  variables  in  a  steady  state  harmonic  system. 

By  far  the  largest  effort  has  been  directed  towards  formulations  which  provide  sym¬ 
metric  coefficient  matrices,  to  again  exploit  the  efficient  algorithms  of  conventional  solver 
systems.  Since  most  of  these  solver  systems  operate  in  simple  arithmetic,  formulations 
which  do  not  produce  a  damping  matrix  are  preferable.  Unfortunately,  the  mixing  of  pri¬ 
mal  (displacement)  and  dual  (pressure)  representations,  which  is  attractive  from  the  point 
of  view  of  minimizing  the  number  of  variables  required  in  the  coupled  problem,  leads  in¬ 
evitably  to  unsymmetric  forms.  Various  methods  of  symmetrizing  such  systems  have  been 
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proposed;  these  can  be  shown  to  be  equivalent  to  reductions  of  three-field  formulations,  in 
which  two  scalar  variables  are  used  in  the  fluid  model,  to  two-field  formulations.  These 
reduction  methods  require  matrix  inversion,  thus  restricting  their  applicability  in  almost  all 
cases.  These  important  limitations  do  not  appear  to  have  been  properly  discussed  in  many 
of  the  formulations  which  have  been  reviewed.  Other  forms  providing  symmetry  result  in  a 
damping  or  damping-like  matrix,  which  presents  its  own  problems. 

The  review  and  development  work  discussed  here  forms  part  of  an  overall  program  of  en¬ 
hancement  of  the  numerical  analysis  capability  within  DREA  for  solution  of  fluid/structure 
interaction  problems.  The  merits  of  a  number  of  the  classical  approaches  are  presented  with 
the  objective  of  establishing  the  formulation (s)  most  compatible  with  the  existing  structural 
finite  element  analysis  system  VAST.47  The  VAST  system  is  fairly  typical  of  conventional 
finite  element  solver  systems;  symmetry  of  the  coefficient  matrices  is  assumed,  and  simple 
arithmetic  is  used  throughout.  Two  well-known  eigenvalue  solution  techniques  are  avail¬ 
able  in  VAST,  the  direct  and  subspace  iteration  algorithms.  Although  these  two  techniques 
are  suitable  in  the  context  of  structural  analysis,  they  are  a  limitation  in  the  application 
of  VAST  to  certain  formulations  of  the  coupled  problem.  It  is  also  not  likely  that  com¬ 
plex  arithmetic  will  be  implemented  in  VAST  in  the  near  future,  hence  those  formulations 
that  result  in  damping  matrices  cannot  currently  be  considered  as  potential  candidates  for 
development. 

Several  of  the  methods  presented  in  this  document  for  treating  coupled  systems  have 
been  evaluated  using  small  test  problems,  and  the  results  are  reported  here.  These  meth¬ 
ods  include  the  structural  analogy  and  displacement-based  approach,  which  require  only 
the  VAST  system,  and  several  of  the  symmetrized  pressure-based  systems,  for  which  a  sep¬ 
arate  program  has  been  developed.  The  latter  is  typical  of  the  reduction  techniques,  being 
relatively  inefficient  and  having  limited  applicability.  The  theoretical  basis  of  this  program 
is  presented  in  the  section  on  scalar  formulations,  and  the  operation  of  the  program  is  pre¬ 
sented  as  an  appendix.  Irrespective  of  its  current  limitations,  the  pressure-based  approach 
embodied  in  this  formulation  can  form  the  basis  of  a  three-field  system  which  should  offer 
a  much  more  efficient  and  general  purpose  capability,  although  a  specialized  eigensolver  is 
required  in  place  of  those  currently  available  in  the  VAST  program.  Development  of  that 
solver  system  is  ongoing;  future  work  will  report  on  its  application  to  practical  problems 
and  performance. 

While  the  pursuit  of  effective  formulations  of  the  undamped  problem  is  worthwhile  in 
other  interaction  problems,  and  makes  for  an  interesting  academic  study  in  the  present 
case,  it  is  well-known  that  damping  plays  a  major  role  in  most  systems  containing  highly 
compressible  fluids.  Damping  can  be  incorporated  into  the  finite  element  formulations  in  a 
number  of  ways,  and  a  section  of  this  report  is  devoted  to  discussion  of  these  approaches. 
Frequency  response  analysis  must  then  generally  be  utilized  to  determine  the  natural  modes 
of  the  coupled  system. 
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In  the  presentation  of  this  work,  it  is  assumed  that  the  reader  is  familiar  with  finite 
element  formulations  of  dynamics  problems.  A  brief  mathematical  development  of  the 
continuum  problem  in  the  individual  fluid  and  structural  domains  is  first  presented  as  a 
basis  for  discussion  of  the  various  numerical  treatments. 


2  Theoretical  Review 


A  mathematical  foundation  for  the  coupled  continuum  problem  is  relatively  easy  to 
develop  and  is  a  necessary  prerequisite  to  a  systematic  consideration  of  the  various  nu¬ 
merical  treatments.  Since  the  formulation  of  the  governing  dynamical  equations  of  the 
structural  system  is  invariant  in  all  of  the  formulations  considered,  this  section  presents  the 
development  of  only  the  fluid  field  equations  in  detail.  These  equations  will  be  referenced 
extensively  in  the  following  s  ons. 

The  treatment  of  the  flu* .  field  embodies  the  classic  assumptions  of  linear  acoustics 
in  which  viscous  and  convective  terms  are  neglected,  and  relations  between  temporal  and 
spatial  distributions  of  pressure,  particle  velocity  and  density  are  sought.  Three  equations 
are  required,  these  typically  comprising  a  continuity  or  conservation  of  mass  relation  in  the 
form  (no  sources), 

0  (1) 

where  p,  is  invariant,  a  momentum  equation, 


P~dt  ~  ~VP  + 

and  a  relation  derived  from  an  adiabatic  equation  of  state 


_  -  OdP 

p  =  p —  = 

Pi 


PtC 

Pi 


in  which  fi  is  the  adiabatic  bulk  modulus.  This  last  equation  can  be  interpreted  as  a  scalar 
constitutive  relation,  since 

dp 

—  =  -V  •  u  (4) 

Pi 


p  =  -pc2V  •  u  (5) 

Hence  stress  in  the  form  of  pressure  is  linearly  related  to  strain  through  a  modulus  fi. 

Appropriate  manipulation  of  these  equations  provides  three  different  functions.  In  terms 
of  density,  we  obtain  a  wave  equation  in  the  form 

^-c2V*p  +  p,V.  g  =  0  (6) 
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Using  equation  (3)  in  the  above,  an  equivalent  pressure  relation  becomes 


(7) 

A  vector  relation  of  some  importance  in  the  current  context  can  be  derived  for  the  equiv¬ 
alent  function  in  terms  of  velocities.  Talcing  the  gradient  and  partial  time  derivative  of 
equations  (1)  and  (2)  respectively,  we  have  in  terms  of  velocity 


-  c2V2p  +  c2p*V  •  g  =  o 


d2v 


dg 


-Qtf  ~  c2 V(V  •  v)  —  =  0 


dt 


or  displacements 


~  c2  V(V  •  u)  -  ^  =  0 
dt2  v  ‘  dt 


(8) 


(9) 


The  necessity  of  assuming  an  irrotational  flow  field  becomes  evident  in  the  above  vector 
relations,  since 

V(V-v)  =  V2v-f  V  *  (V  *  v)  (10) 


Irrotationality  of  the  flow  (or  displacement)  is  required  to  eliminate  the  cross-product  term. 
With  such  an  assumption,  equations  (8)  and  (9)  can  be  cast  in  scalar  form  via  the  intro¬ 
duction  of  velocity  and  displacement  potential  functions  respectively. 

It  can  be  shown  that  the  contributions  of  gravity  effects,  such  as  surface  waves,  are 
generally  small  in  comparison  to  compressibility  contributions  in  a  coupled  system  and 
the  gravitational  terms  in  the  above  relations  are  often  neglected.  This  term  is  obviously 
retained  in  the  formulation  of  solutions  of  the  sloshing  problem,  in  static  problems  and  in 
problems  of  large  scale  in  which  variation  of  pre-stresses  may  be  important.  The  common 
forms  of  the  scalar  equations  are  the  classical  wave  equations  of  potential  or  pressure  in  the 
fluid  domain 


V2<f> 


aV 

:J  dt 2 


(11) 


V2p 


a2p 

tdt2 


(12) 


In  the  remainder  of  this  report,  the  symbol  p  will  replace  p ,,  except  where  the  distinction 
is  required. 

The  boundary  conditions  for  the  fluid  field  are  reasonably  self-evident.  For  the  velocity 
potential  and  pressure  formulations  respectively,  they  are  appropriate  combinations  of 


d<f>  dp 

^  =  u''n  dii  =  ~pUt ' n 


(13) 
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on  the  elastic  interface, 


(14) 


on  a  rigid  wall,  and 


p4>  =  o 


P  =  0 


(15) 


on  a  free  surface  of  zero  ambient  pressure.  Here,  u,  is  the  structural  displacement  vector, 
and  the  overdot  notation  is  used  for  convenience  to  represent  time  differentiation. 

An  equivalent  scalar  formulation  using  a  displacement  potential  ip  in  place  of  the  velocity 
potential  can  also  be  developed,  the  appropriate  boundary  conditions  becoming 


Vip  ■  n  =  u,  •  n 

(16) 

on  the  elastic  interface, 

Vip  •  n  =  0 

(17) 

on  a  rigid  wall,  and 

ip  =  0 

(18) 

on  a  free  surface.  Other  impedance  type  boundary  conditions  are  sometimes  required,  such 
as  the  radiation  condition  for  unbounded  domains,  and  treatments  of  partially  transmitting 
boundaries.  These  will  be  discussed  in  more  detail  in  following  sections. 

The  preceding  development  provides  the  differential  equations  and  boundary  conditions 
relevant  to  the  complete  fluid  domain.  To  utilize  the  finite  element  method,  an  integral 
relation  is  required,  which  when  solved,  provides  an  approximate  solution  to  these  equations 
and  boundary  conditions.  Two  related  procedures  are  used  to  obtain  the  integral  functional: 
direct  application  of  the  weighted  residual  method,  and  formulation  of  the  system  on  the 
basis  of  a  variational  principle.  The  former  is  much  more  general,  but  there  is  always 
impetus  to  establish  a  variational  principle,  if  one  exists,  since  it  imparts  attractive  features 
such  as  bounding  to  a  well-posed  numerical  solution.  In  this  case,  variational  principles 
based  on  energy  can  be  developed  for  both  the  fluid  and  structural  domain,  provided  that 
the  individual  systems  considered  are  conservative.  For  such  cases,  the  minimization  of  the 
functional  corresponds  exactly  to  the  solution  of  the  relevant  differential  equations  with 
their  proper  boundary  conditions.  In  general,  if  a  variational  principle  exists,  it  can  often 
be  obtained  using  appropriate  operations  on  the  basic  weighted  residual  statement,  as  is 
presented  below.  Here,  as  is  often  most  logical,  the  Galerkin  weighted  residual  method  is 
utilized,  in  which  the  weighting  functions  and  interpolation  functions  are  the  same. 

For  the  scalar(p,^)  relations  developed,  the  general  form  of  the  basic  weighted  residual 
statement  is 
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f  6p(Vp  +  pii,)  •  n)  dT  =  0 
J  r. 


(19) 


where  ft/  represents  the  entire  fluid  volume,  and  T*,  represents  the  flexible  portion  of  the 
fluid/structure  interface.  Assuming  that  the  variation  of  the  pressure  function  6p  satisfies 
the  boundary  condition  of  equation  (15),  (the  essential  boundary  condition),  equation  (19) 
can  be  integrated  once  by  parts  to  obtain  a  weak  form  of  the  governing  equations  in  the 
form  of  a  variational  statement 

6  f  ^(pJ  4-  pc2Vp  •  Vp)  dft  -  f  ppii ,  •  n  dT  =0  (20) 

Jnf  2  J r„ 

The  volume  integral  in  this  functional  can  be  used  as  the  basis  for  finite  element  solutions 
of  rigid-walled  acoustic  cavity  problems.  In  that  case,  the  absence  of  a  surface  integral 
imposes  the  boundary  condition  of  equation  (14)  automatically.  Here,  the  surface  integral 
enforces  the  coupling  boundary  condition  at  the  elastic  interface.  Equation  (20)  is  in  fact  a 
variational  principle,  in  which  the  kernel  of  the  volume  integral  is  the  acoustic  Lagrangian. 
The  surface  integral  accounts  for  the  work  done  by  the  external  loads  as  seen  by  the  fluid. 

While  the  formulation  of  the  scalar  relations  in  terms  of  a  velocity  potential  is  prob¬ 
ably  more  familiar,  the  use  of  a  displacement  potential  is  usually  more  convenient  in  the 
present  application,  since  a  discrete  system  can  be  obtained  without  the  presence  of  velocity 
dependent  terms.  This  point  will  be  discussed  further  in  the  next  section. 

Following  the  same  process  as  described  in  some  detail  for  the  fluid,  a  weighted  residual 
process  can  also  be  used  to  develop  a  variational  form  of  the  governing  equations  of  the 
structural  system.  In  this  case,  the  functional  contains  terms  representing  the  kinetic  and 
potential  energies,  the  latter  comprising  the  strain  energy,  work  of  the  body  forces,  and 
work  done  by  external  loads.  The  fluid  pressure  term  is  separated  from  other  loads,  and 
is  the  source  of  the  fluid/structure  coupling.  A  typical  variational  form  in  terms  of  the 
structural  displacements  is 

&  f  (i/>,u*u,  -  ^VuJ'DVu, -Xru,)  dft-  f  Fru ,  dT  -  f  (u, -n)pdT  =  0  (21) 

Jo,  22  j p,  Jrm 

where  D  is  the  elasticity  matrix,  X  are  the  body  forces,  and  F  are  the  external  loads.  The 
solution  of  the  coupled  problem  thus  involves  the  simultaneous  satisfaction  of  relations  in 
the  form  of  (20)  and  (21).  The  finite  element  solution  proceeds  in  the  standard  manner 
of  assuming  local  approximating  functions  for  the  (u ,,p,<p,ip)  variables,  as  appropriate, 
numerically  evaluating  the  integrals  and  assembling  the  elemental  matrices  into  an  overall 
system.  The  following  section  discusses  the  numerous  types  of  discrete  formulations  of  these 
relations  which  have  been  proposed. 

Within  the  assumptions  of  linear  acoustics  and  small  displacement  theory,  there  is  no 
mathematical  approximation  necessary  in  the  derivation  of  the  above  functionals.  Of  course, 
any  numerical  procedure  can  introduce  errors  into  the  solution  and  destroy  the  features  of 
such  variational-based  solutions.  For  example,  the  use  of  reduced  integration,  which  has 
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been  found  beneficial  in  certain  formulations,  may  corrupt  the  underlying  functional,  and 
lead  to  a  solution  of  an  unknown  set  of  governing  equations. 

3  Discrete  Formulations  of  the  Coupled  Eigenproblem 

The  numerous  discrete  formulations  of  the  coupled  problem  can  be  crudely  divided  into 
three  groups  based  on  the  main  objective  of  the  approach.  One  group  is  concerned  pri¬ 
marily  with  the  response  of  the  structure,  another  is  concerned  primarily  with  the  acoustic 
response  of  the  fluid,  and  a  third  group  treats  both  domains.  In  the  case  of  the  former 
two  groups,  simplifying  assumptions  about  the  behaviour  of  the  complementary  domain 
can  often  be  exploited  in  order  to  substantially  decouple  the  problem.  This  section  begins 
with  a  very  brief  discussion  of  these  first  two  types  of  approaches,  but  the  emphasis  will  be 
on  formulations  which  consider  the  response  of  both  structured  and  fluid  domains.  Further 
discussion  of  uncoupled  analysis  techniques  can  be  found  in  several  of  the  references  quoted. 

Much  of  the  early  finite  element  work  in  the  area  of  fluid/structure  coupling  was  directed 
toward  the  prediction  of  the  effect  of  the  fluid  on  the  dynamic  response  of  large  structures 
such  as  gravity  dams.2,38,45  In  these  cases,  only  the  structural  response  in  the  low  frequency 
range  (below,  say,  20  Hz)  is  of  interest.  For  these  problems,  ( ka  small),  the  fluid  effect  is 
almost  purely  reactive,  and  compressibility  effects  could  in  most  cases  be  neglected.  This 
leads  to  the  well-known  added  mass  concept,  and  if  free  surface  gravity  waves  are  included, 
to  an  added  damping.  Finite  elements  of  various  forms  have  been  applied  with  much  success 
to  this  problem,  and  although  it  is  essentially  an  unbounded  one,  these  developments  have 
laid  the  groundwork  for  many  of  the  finite  element  formulations  of  the  bounded  problem. 
Reference  45  remains  essential  reading  for  anyone  working  in  this  area. 

Similar  developments  occurred  during  roughly  the  same  time  period  in  the  use  of  finite 
elements  for  acoustic  modal  analysis  and  the  prediction  of  acoustic  response  of  coupled 
systems.4,6,14,37  To  a  large  extent,  the  structural  response  was  of  secondary  importance  in 
these  studies,  and  approximate  structural  modes  were  often  assumed  representative  of  the 
boundary  forcing  functions  in  the  acoustic  domain.  Here  too,  finite  element  formulations, 
usually  pressure-based,  have  proven  very  effective  for  the  prediction  of  acoustic  response 
at  low  to  moderate  frequencies.  The  modal  approximation  approaches  have  a  number  of 
merits,  and  can  be  used  effectively  for  complete  solution  of  the  coupled  problem,  as  will  be 
discussed. 

Numerical  treatments  which  utilized  finite  element  modelling  of  both  domains  began 
to  appear  in  the  mid-1970’s,  although  the  mathematical  basis  for  their  development  had 
probably  been  presented  a  decade  previously  in  the  work  of  Gladwell  and  Zimmerman,30,21 
and  again  by  Zienkiewicz.46  Since  then,  a  continuous  discussion  of  this  problem  can  be 
found  in  the  literature. 
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It  is  instructive  to  divide  the  various  finite  element  formulations  into  subsections  based 
on  the  particular  variable  or  variables  used  to  characterize  the  fluid  field.  In  this  way, 
the  merits  and  drawbacks  of  each  general  method  can  be  discussed,  and  the  links  to  other 
approaches  made  in  some  consistent  fashion.  A  brief  summary  of  the  relevant  characteristics 
is  also  given  at  the  end  of  the  section.  A  pressure  formulation  is  presented  in  some  detail 
to  define  the  basic  terms  that  occur  in  some  form  in  almost  all  of  the  approaches. 

The  choice  of  field  variables  can  have  different  effects  on  solutions  of  transient  and  steady 
state  problems.  In  this  discussion,  we  consider  only  the  ramifications  for  the  solution  of 
the  eigenproblem.  The  application  of  some  of  the  formulations  as  applied  to  time  domain 
solutions  cam  be  found  in  certain  of  the  references. 1'4*16,19,22,41 

3.1  Scalar  Formulations 

The  use  of  scalar  variables  to  describe  the  fluid  field  is  obviously  appealing  since  it 
results  in  a  minimum  number  of  fluid  degrees  of  freedom,  while  automatically  satisfying 
the  irrotationality  condition  arising  from  the  small  displacement  assumption  in  the  fluid. 
Thus,  these  formulations  are  immediately  applicable  for  solution  of  the  scalar  wave  or 
Helmholtz  equation.  Numerous  variations  have  been  popular,  most  notably  those  utilizing 
fluid  pressure,  velocity  potential  and  more  recently  displacement  potential.  The  scalar 
analogy  is,  as  its  title  suggests,  analogous  to  these  formulations,  but  it  is  discussed  separately 
here  because  of  its  rather  unique  basis  and  requirements. 

3.1.1  Pressure 

The  description  of  the  fluid  field  in  terms  of  a  pressure  variable  is  conceptually  appealing 
as  well  as  practical  from  the  viewpoint  of  applying  certain  boundary  conditions.  This  ap¬ 
proach  can  be  used  very  effectively  in  cavity  acoustics  problems  and  uncoupled  systems  such 
as  the  prediction  of  added  mass,  but  poses  some  difficulties  for  compressible  fluid/structure 
interaction. 

To  start,  the  coupled  equations  arising  from  the  variational  statements,  equations  (20) 
and  (21),  or  more  generally  from  the  application  of  a  Galerkin  finite  element  approach  if 
damping  is  to  be  included,  can  be  cast  in  the  form  of  two  discrete  differential  equations.  In 
keeping  with  standard  structural  finite  element  methodology,  all  of  the  approaches  reviewed 
utilize  the  set  of  nodal  structural  displacements  6  as  one  set  of  unknowns.  In  terms  of  the 
set  of  nodal  fluid  pressures  p,  these  two  equations  are 

lM}6  +  [D.)6  +  lK)8  =  [L)p+{Ft} 

\Q}p+[Df]p  +  \H]p=  -p[L\T6 
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(22) 

(23) 


with  p  the  fluid  density,  and  Ft  the  vector  of  structural  forces  independent  of  the  fluid. 
The  structural  equation  contains  the  familiar  mass,  damping  and  stiffness  matrices  for  a 
multi-degree  of  freedom  system.  The  coupling  terms  occur  as  forcing  functions  in  each 
equation;  for  the  fluid,  the  acceleration  of  the  structural  interface  produces  fluid  pressures, 
and  for  the  structure,  those  induced  pressures  produce  structural  forces.  In  the  standard 
approach,  the  fluid  and  structure  have  been  discretized  into  finite  elements  within  which 
the  structural  displacements  and  fluid  variables,  here  pressure,  are  described  by 

*  =  (24) 

P  =  T,NiPi  (25) 

and  the  IV,-  represent  interpolation  functions,  which  can  be  different  between  the  fluid  and 
the  structure.  The  subscripted  field  variables  are  the  nodal  values,  and  the  superscripts 
denote  the  variable  with  which  the  interpolation  function  is  associated.  With  reference 
to  equation  (20),  the  fluid  coefficient  matrices  in  equation  (2)  are  assembled  forms  of  the 
inertial  terms 

«  =  ?  /„,  N’N’ dn + ;  /„,  dn  <26> 

the  fluid  damping  terms 

(<*/>« = v  L  N’N’ da  <27> 

c  Jnf  j 

and  the  acoustic  ‘stiffness’  terms 

htj  =  [  VATf  •  VNf  d(i  (28) 

Jn,  3 

The  L  matrix  is  the  coefficient  matrix  relating  pressure,  acting  normal  to  the  interface,  to 
global  forces  at  the  nodes 

(29> 

Neglecting  the  damping  terms  for  the  moment  (radiation  damping  is  usually  assumed  small, 
but  concentrated  acoustic  damping  may  be  quite  significant),  the  eigenvalue  problem  asso¬ 
ciated  with  equations  (l)  and  (2)  becomes 


Despite  the  inherent  symmetry  in  the  K,  H  and  Q  submatrices,  the  assembled  form  of  this 
coupled  system  of  this  formulation  is  unsymmetric.  Conventional  solution  systems  such  as 
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those  available  in  the  VAST  program  cannot  be  used  to  solve  this  problem.  Equation  (30) 
is  the  basis  of  a  host  of  variations  on  the  theme  of  scalar  formulations. 

The  lack  of  symmetry  in  this  system  is  a  direct  result  of  the  mixed  nature  of  this  common 
formulation,  which  utilizes  displacement  variables  in  the  structure  and  dual  (force)  variables 
in  the  fluid.  On  a  mathematical  basis,  the  finite  element  models  in  the  two  domains  are 
approximating  solutions  of  quite  different  quantities;  a  displacment-based  solution  generally 
satisfies  compatibility,  but  approximates  equilibrium,  while  complementary  formulations 
satisfy  equilibrium,  but  approximate  compatibility.  Thus,  as  pointed  out  by  Tabarrok,43 
the  only  consistent  solution  is  the  exact  solution,  which  is  never  obtained.  The  consequences 
of  such  mixed  formulations  are  often  difficult  to  assess  a  priori  ,  but  at  the  very  least,  the 
solutions  cannot  be  used  as  bounds,  as  is  often  the  case  for  consistent  formulations. 

The  inconsistency  can  be  avoided  if  a  complementary  formulation  utilizing  force  vari¬ 
ables  is  used  for  the  structure.  To  this  author’s  knowledge,  no  results  of  finite  element 
formulations  utilizing  the  consistent  dual  principle  have  appeared  in  the  literature. 

Because  the  undamped  coupled  problem  can  be  developed  from  a  variational  principle 
in  quadratic  form,  symmetric  forms  of  equation  (30)  must  exist.  In  fact,  four  independent 
modifications  of  that  system  can  be  made  in  order  to  symmetrize  the  system,  each  of  which 
involves  a  substitution  and  requires  an  inverse  of  one  of  the  four  submatrices  K,M,H  or  Q. 
For  example,  a  useful  modified  relation  is  obtained  if  the  substitution 

p  =  H~l  {pLT  S  +  Qp)  (31) 


is  made  in  the  structural  partition  of  equation  (30).  Multiplication  of  equation  (31)  by  Q 
and  division  by  p  gives  a  new  second  row,  and  a  standardized  system  is  obtained  in  the 
form 

M  +  pLH~lLr  LH-'Q 

IQH-'Q 


T  K  o' 

2 

.1  0  }<3  , 

-  OJ 

]]{;}-{:) 


(32) 


where  a  dash  implies  symmetry  across  the  diagonal.  Application  of  the  three  alternative 
substitutions  which  can  be  used  to  obtain  symmetric  forms  of  equation  (30)  provide  the 
systems 


T  K  +  pLQ 

~lLT 

IL 

rr  m  o 

2 

ii° 

-  U) 

-LQ~XH 


MK~lM 
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M 

0 


°  iim  _/ o\ 
iH\\  l  p  J  1  0  J 


-MK~lL 
\Q  +  LTK~lL 


1  /  6  \  _  /  0  1 
\  P  J  \  0  J 


(33) 


(34) 
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(35) 


KM~lK 


-KM~XL 
\H  +  LrM~1L 


or 


K  0 
0  )Q 


1  /  8  1  _  /  0  1 
W  p  j  \  °  j 


Each  of  the  equations  (32)  to  (35)  is  a  symmetric  system  which  may  be  soluable  with 
conventional  finite  element  solution  techniques,  provided  that  the  system  can  in  fact  be 
formulated.  Since  several  of  these  systems  have  received  attention  in  the  literature,  it  is 
worthwhile  to  consider  the  requirements  and  limitations  in  each  approach,  a  point  not  well- 
considered  in  publications  reviewed  by  this  author.  Taking  each  equation  in  turn,  we  note 
that  system  (32)  requires  a  positive  definite  structural  stiffness  matrix  K  and  fluid  matrices 
H  and  Q.  In  most  cases,  the  rigid  body  constraint  of  the  structure  is  not  problematic,  and 
Q  is  non-singular;  however,  for  closed  cavities,  the  H  matrix  iB  singular,  and  inversion  is 
possible  only  through  the  constraint  of  at  least  one  pressure  degree  of  freedom.  For  such 
problems,  experience  has  shown  that  such  a  constraint  always  corrupts  the  solution.  More 
will  be  said  shortly  about  this  particular  system. 

System  (33)  requires  the  inversion  of  the  fluid  mass  matrix  Q,  which  is  always  realizable, 
and  which  will  not  be  particularity  onerous  if  a  diagonalized  matrix  can  be  utilized,  or  if  the 
fluid  domain  is  small.  Unfortunately,  the  presence  of  a  singular  H  matrix  in  the  diagonal 
block  of  the  assembled  stiffness  matrix  renders  that  matrix  singular.  Shifting  will  not 
alleviate  this  problem,  since  the  shift  involves  the  use  of  the  singular  H  matrix.  Thus, 
equation  (33)  is  not  a  useable  formulation  for  the  general  problem. 

System  (34)  is  attractive  in  that  it  retains  the  banded  nature  of  the  individual  systems  in 
the  assembled  stiffness  matrix,  and  that  the  use  of  shifting  may  allow  the  presence  of  singular 
H  matrices.  The  approach  requires  a  positive  definite  mass  matrix,  and  the  inversion  of  the 
structural  stiffness  matrix.  Without  considering  the  former,  the  approach  must  be  dismissed 
on  the  premise  that  the  such  an  inversion  is  not  practical  for  anything  other  than  purely 
academic  problems.  This  approach  has  apparently  been  successfully  utilized  in  References 
10  and  11  for  some  very  small  systems. 

The  final  system  of  equation  (35)  requires  a  structural  mass  matrix  inversion,  which  is 
tractable  for  diagonalized  forms.  As  well,  the  application  of  shifting  should  allow  solutions 
of  singular  fluid  systems.  Aside  from  the  poorly  banded  form  of  the  assembled  stiffness 
matrix,  this  approach  presents  several  other  numerical  difficulties.  First,  the  structural 
mass  matrix  may  not  be  well-conditioned  when  elements  utilizing  rotational  degrees  of 
freedom  are  used  in  the  structural  model.  This  characteristic  arises  from  the  disparate 
values  of  translational  and  rotational  mass  components,  and  is  not  alleviated  through  the 
natural  resort  to  a  diagonal  mass  matrix.  Since  these  element  types  are  a  necessary  part  of  a 
realistic  structural  model,  this  is  a  non-trivial  problem.  Second,  the  squaring  of  terms  in  the 
structural  stiffness  matrix  can  lead  to  numerical  overflows  when  structural  constraints  are 
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imposed  through  the  use  of  luge  spring  stiffnesses,  such  as  are  utilized  in  current  versions  of 
VAST.  The  use  of  smaller  stiffness  values  may  adversely  affect  the  accuracy  of  the  prediction 
of  structural  response.  As  well,  there  is  still  considerable  matrix  manipulation  required  even 
when  a  diagonal  mass  matrix  is  utilized.  Presumably,  the  removal  of  the  constrained  degrees 
of  freedom  would  eliminate  this  second  problem,  but  the  first  would  remain.  Further  work  is 
required  to  investigate  if  rotational  masses  could  be  modified  to  provide  better  conditioning 
of  the  structural  mass  matrix. 

Despite  its  inherent  limitations,  the  reduced  formulation  in  equation  (32)  is  straightfor¬ 
ward  to  produce,  and  can  be  used  in  the  solution  of  constrained  types  of  coupled  problems. 
This  system  has  been  coded  in  a  stand-alone  computer  program  COUPLE,  which  creates 
the  fluid  submatrices  Q,H ,  and  L,  carries  out  the  required  manipulations,  and  assembles 
the  overall  system  with  the  structural  submatrices  K  and  M.  This  program  incorporates 
the  previously  developed  acoustic  finite  element  generator  program  ACOUSTIC.44  The 
program  is  also  configured  to  construct  the  system  of  equation  (35),  although  the  practi¬ 
cal  application  of  this  approach  is  still  quite  limited  due  to  the  problems  discussed  above. 
Further  discussion  of  this  program  is  contained  in  the  appendix. 

Of  note  in  equation  (32)  is  the  conventional  added  mass  effect  provided  by  this  sym¬ 
metrized  formulation  in  the  limit  of  fluid  incompressibility.  As  the  sound  speed  tends  to 
infinity,  all  terms  involving  Q  vanish  (aside  from  the  small  surface  wave  contributions,  see 
equation  (26)),  and  equation  (32)  reduces  to 

[[*]  -w2[M  +  MA]]{S}  =  0  (36) 

in  which 

\MA\  =  pLH~lLT  (37) 

This  is  identical  to  the  form  utilized  in  the  current  VAST  fluid  element  formulation  for 
added  mass, 23  and  shows  how  the  limiting  assumption  of  incompressibility  uncouples  the 
general  system. 

As  first  pointed  out  by  Morand  and  Oh  ay  on30,  the  symmetrized  systems  of  equa¬ 
tions  (32)  to  (35)  can  be  interpreted  as  statically  reduced  forms  of  a  three-field  system. 
For  example,  equation  (32)  is  the  reduced  form  of  an  expanded  symmetric  system  in  the 
(6,p,ip)  variables,  in  which  ip  can  be  shown  to  be  equivalent  to  a  displacement  potential 
function  in  the  fluid.  Three-field  systems  offer  a  number  of  advantages  over  their  two-field 
counterparts,  and  will  be  discussed  further  in  the  subsection  covering  that  topic.  The  relar 
tionships  between  the  two  and  three-field  systems  are  discussed  in  some  detail  by  Fellipa, 17,18 
who  also  considers  at  least  to  some  degree  the  utility  of  the  various  reduced  systems,  which, 
as  pointed  out,  have  quite  serious  limitations. 
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3.1.2  Displacement  Potential 


A  displacement  potential  formulation  for  the  fluid  domain  can  be  developed  from  a  direct 
application  of  the  Galerkin  weighted  residual  technique  to  the  wave  equation  in  terms  of 
displacement  potential.  A  weak  form  of  the  system  becomes 


+  c2vv>  •  v^)  dn  -  c2 


with  a  discrete  system  obtained  from  a  variation  with  respect  to  tp, 


[. Hm+[Qm+[L]T{6 >  =o 


(38) 


(39) 


The  coupling  term  in  the  otherwise  unchanged  structural  equation  can  be  obtained 
through  substitution  of  the  constitutive  relation  of  equation  (3),  and  the  displacement 
potential  wave  equation  into  the  last  term  of  equation  (21), 


-  I  (u,  •  n)p  dT  ~  p  f  (u,  •  n)tp  dT 
Jrm  J r» 

The  assembled  system  then  becomes 

[[?;]-■[; «]]{:)-(:) 


(40) 


(41) 


Evidently,  the  use  of  a  displacement  potential  does  not  result  in  a  symmetric  formulation. 
The  form  of  the  assembled  system  is  very  similar  to  the  pressure  formulation  of  equa¬ 
tion  (30),  and  there  is  no  apparent  advantage  in  using  such  an  approach  on  its  own.  The 
reversal  of  the  coupling  terms  suggests  that  combinations  of  the  two  variables  might  yield 
some  form  of  symmetric  system.  As  was  alluded  to  earlier,  the  utility  of  the  displacement 
potential  variable  indeed  becomes  evident  in  the  three-field  formulations. 


3.1.3  Velocity  Potential 

A  reformulation  of  the  basic  scalar  problem  using  a  fluid  velocity  potential  in  place  of 
the  fluid  pressure  has  been  suggested  by  Everstine.15  The  relationship  between  the  pressure 
and  velocity  potential  is  that  commonly  associated  with  potential  flow  theory 

P  «  (42) 

Making  this  substitution  in  equation  (30),  dividing  the  fluid  partition  of  that  equation  by 
-p  and  integrating  in  time  gives  the  system 
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Regardless  of  the  presence  of  fluid  or  structural  damping  terms,  the  time  integration  required 
in  this  formulation  moves  the  coupling  terms  in  L  to  the  damping  matrix.  For  the  eigenvalue 
problem,  this  transforms  the  undamped,  unsymmetric  problem  into  a  complex  symmetric 
one.  For  systems  which  operate  in  simple  arithmetic,  this  is  not  a  viable  formulation.  For 
other  types  of  time  domain  analyses,  the  approach  does  have  redeeming  features. 


3.1.4  Scalar  Analogy 

An  alternative  method  of  generating  the  required  fluid  submatrices  H  and  Q  in  the 
coupled  problem  is  to  use  standard  structural  solid  elements  with  three  displacement  degrees 
of  freedom  per  node,  in  what  Everstine  has  termed  a  fluid-structural  analogy.16  In  this  quite 
elegant  trick,  the  wave  equation 

vv=4^  (44) 

can  be  solved  with  structural  formulations  based  on  the  Navier  equations  of  elasticity  in 
the  form 

(A  )  V2Ul  +  (“2 ,*v  +  “3.")  =  (45) 

when  the  uj  and  uj  components  of  displacement  are  zero,  and  the  conditions 

X  +  2M  _  t 


P_l_ 

H  c 2 

are  satisfied.  Everstine  has  related  these  requirements  for  the  Laml  constants  of  elasticity, 
H  (the  shear  modulus)  and  A,  to  the  more  conventional  material  property  constants  E  and  u, 
which  are  more  easily  specified  in  finite  element  programs  as  input  data.  A  typical  approach 
is  to  arbitrarily  assign  n  a  value  of  unity  and  specify  Young’s  modulus  E  large,  say  1.0E+15. 
The  specification  of  a  Poisson’s  ratio  v—E/2  and  a  density  of  1/c3  then  makes  the  solution 
of  equations  (45)  and  (44)  equivalent.  An  uncoupled  acoustic  cavity  analysis  capability  is 
then  immediately  available  from  an  eigenvalue  analysis  using  stiffness  and  mass  matrices 
assembled  from  standard  elements.  An  added  fluid  mass  prediction  capability  is  also  easily 
defined  by  using  standard  constraints  on  the  exterior  surface  of  a  fluid  model. 

The  advantage  of  this  approach  is  that  no  special  elements  need  be  developed,  and 
standard  pre/post  processors  can  be  used  for  interpretation  of  input  and  output  data; 
however,  for  a  program  such  as  VAST,  there  are  several  distinct  disadvantages.  First,  two 
degrees  of  freedom  at  each  fluid  node  must  be  constrained  to  simulate  a  scalar  field.  These 
constrained  freedoms  are  not  eliminated  in  the  VAST  solution,  hence  the  problem  size  can 
be  significantly  larger  than  necessary. 
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While  quite  elegant  for  uncoupled  analysis,  the  major  drawback  in  the  present  context 
is  that  the  scalar  formulation  of  the  coupled  problem  is  in  exactly  the  unsymmetric  form 
given  previously  for  the  pressure  formulation.  Since  several  well  developed  capabilities 
already  exist  for  added  mass  prediction  and  the  programs  ACOUSTIC  or  COUPLE  are 
available  for  efficiently  generating  the  matrices  required  for  acoustic  cavity  modal  analysis, 
the  scalar  analogy  does  not  offer  any  enhancement  to  the  current  capability.  Finally,  the 
Poisson’s  ratio  specification  required  in  this  analogy  is  structurally  unrealistic,  and  causes 
problems  in  programs  such  as  VAST  which  check  for  such  characteristics.  The  analogy  has 
been  tested  with  the  VAST  program,  and  a  recent  comparison  with  results  from  acoustic 
elements  is  given  in  the  next  section.  Further  discussion  of  this  approach  in  the  context  of 
VAST  operation  can  also  be  found  in  Reference  3. 

3.2  Vector  Formulations  -  Displacement  Fluid  Variables 

Since  virtually  all  structural  analysis  codes  use  displacements  as  primary  variables,  it 
is  natural  to  formulate  the  coupled  problem  completely  in  terms  of  displacements.  Such 
an  approach  should  avoid  the  inconsistency  inherent  in  mixed  formulations  and  thus  lead 
directly  to  symmetric  forms.  This  is  in  fact  the  case,  although  the  immediately  appar¬ 
ent  drawback  is  the  larger  number  of  fluid  variables  required.  The  disadvantage  of  larger 
problem  size  may  be  offset  by  the  bandedness  which  can  usually  be  conserved  in  consistent 
formulations. 

A  less  obvious  but  serious  drawback  is  that  a  straightforward  application  of  the  displace¬ 
ment  variable  formulation  will  not  satisfy  the  irrotationality  condition  which  is  implicit  in 
the  use  of  the  scalar  formulations.  The  result  is  the  presence  of  what  have  been  called 
circulation  modes  in  the  fluid  domain.  These  modes  are  displacement  configurations  which 
involve  zero  volumetric  strain  in  the  elements,  their  exact  form  and  frequency  often  being 
dependent  on  the  particular  grid,  element  order  and  integration  point  location.  In  general, 
these  modes  will  have  non-zero  frequencies,  and  for  detailed  grids,  the  number  of  circulation 
modes  may  make  determination  of  the  true  compressible  modes  very  difficult. 1,24,27 

Although  a  simple  solution  would  appear  to  be  the  use  of  the  form 

^-c’V’u  =  0  (48) 

as  a  starting  point,  a  weak  form  of  this  system  produces  untractable  surface  integrals. 
Recalling  that  this  equation  is  a  reduced  form  of  equation  (9),  in  which  the  cross-product 
term  of  the  expansion  (10)  was  assumed  to  vanish,  a  more  general  starting  point  is 


Unfortunately,  there  is  now  no  guarantee  that  the  displacement  field  satisfies  the  irrotar 
tionality  requirement.  This  is  the  form  that  an  unmodified  finite  element  program  such 
as  VAST  will  emulate.  One  solution  to  this  problem  is  the  addition  of  a  penalty  function 
to  equation  (49)  based  on  the  relation  of  equation  (10).  In  this  case,  the  starting  point 
becomes 

—  eJ  V(V  •  u)  +  qcj V  *  (V  *  u)  =  0  (50) 

with  a  a  non-zero  penalty  coefficient.  A  weak  form  of  this  equation  can  be  obtained  (two 
rules  of  vector  differentiation  are  required)  in  the  form 


+  (V  •  $u)(V  •  u)  +  ac2(V  *  $u)(V  *  u)  )  dfl  -I-  j  *u(V  •  u) 


dT  =  0 


(51) 


The  surface  integral  produces  a  load  term  in  both  domains  which  is  a  function  of  the 
strain  at  the  interface.  For  rigid-walled  cavities,  and  coupled  systems  in  which  consistent 
interpolation  functions  are  utilized  in  the  two  domains,  the  integral  vanishes.  If  mixed 
interpolants  are  used,  a  fictitious  load  term  results. 

When  utilized,  the  discrete  form  of  the  penalty  function  takes  the  form  of  an  added 
stiffness  term  Hr.  The  assembled  form  of  the  coupled  system  becomes 
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The  use  of  large  values  of  a,  on  the  order  of  one  to  ten  times  the  fluid  bulk  modulus, 
suppresses  the  circulation  modes  in  most  cases,  and  does  not  appear  to  adversely  degrade  the 
prediction  of  the  true  modes.1’24  In  current  versions  of  the  program  ADINA48,  which  offers 
displacement-based  formulations  with  or  without  rotational  stiffness  terms,  the  penalty 
coefficient  is  taken  as  one  thousand  times  /?. 

Certain  displacement  formulations19  also  include  interface  elements  between  the  fluid 
and  structure  in  order  to  decouple  the  tangential  displacements,  consistent  with  the  inviscid 
assumption.  This  decoupling  has  been  shown  to  provide  better  solutions  than  a  completely 
coupled  form  in  some  cases,  which  appears  to  be  too  rigid.  Another  suggestion  for  tangential 
decoupling  in  cases  where  interface  elements  are  not  available,  is  the  use  of  short  stiff  rod 
elements  between  fluid  and  structural  nodes  at  the  interface.1 

To  utilize  this  formulation  with  a  conventional  displacement-based  finite  element  pro¬ 
gram  such  as  VAST,  a  modified  elasticity  matrix  is  required  for  the  fluid  elements  in  order 
to  enforce  the  zero  shear  stress  model.  For  3-dimensional  problems,  the  standard  D  matrix 
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must  take  the  form 


'  p  p  p  0  0  0  ' 

p  P  P  0  0  0 

P  0  P  0  0  0  .  . 

oooooo 

0  0  0  0  0  0 

.0  0  0  0  0  0. 

where  P  is  the  bulk  modulus  of  the  fluid.  Appropriately  reduced  forms  of  this  matrix  are 
utilized  for  2-dimensional  or  axisymmetric  analyses.  The  VAST  program  currently  allows 
direct  input  of  the  elasticity  matrix  for  certain  membrane,  plate,  shell  and  solid  elements. 

The  consistency  of  the  displacement  variable  approach  is  appealing,  since  any  graphics 
and  post-processing  features  developed  for  the  general  analysis  program  can  readily  be  uti¬ 
lized.  The  method  has  been  applied  in  the  analysis  of  several  simple  problems,  the  results 
of  which  are  presented  in  the  next  section.  The  formulation  can  be  used  to  solve  the  cou¬ 
pled  problem  in  which  no  constraints  are  applied  in  the  fluid,  for  which  the  symmetrized 
scalar  approaches  were  unworkable  due  to  the  presence  of  the  singular  fluid  stiffness  ma¬ 
trix.  In  both  coupled  and  uncoupled  solutions,  circulation  modes  are  evident,  and  it  is 
suggested  that  these  would  be  very  problematic  in  large  scale  analysis.  The  inclusion  of  a 
penalty  stiffness  term  for  selected  elements  would  not  be  a  difficult  modification  in  VAST, 
and  the  existing  interface  elements  could  readily  be  used  to  provide  tangential  decoupling. 
Additional  capability  for  entering  the  elasticity  matrix  would  be  required;  again  this  would 
require  only  minor  modification. 

The  application  of  this  method  in  the  solution  of  general  problems  is  not  without  dif¬ 
ficulties,  and  the  simple  configurations  considered  in  the  next  section  highlight  some  of 
these.  First,  the  highly  singular  nature  of  the  elasticity  matrix  imposes  stringent  require¬ 
ments  on  the  boundary  conditions  which  must  be  applied  to  the  problem  in  order  to  avoid 
ill-conditioning  of  the  overall  stiffness  matrix.  For  general  geometries,  the  use  of  awkward 
skewed  coordinates  is  also  necessary  in  the  application  of  correct  boundary  conditions.  As 
well,  even  when  a  rotational  stiffness  is  available,  it  is  difficult  to  specify  an  appropriate 
value  of  the  penalty  coefficient  a  priori,  and  repeated  evaluations  and  plotting  of  mode 
shapes  may  be  required  to  determine  true  modes.  In  Reference  35,  poor  performance  of 
the  displacement-based  method  was  also  noted  in  the  prediction  of  coupled  eigenmodes  for 
nearly-incompressible  fluids.  These  immediately  apparent  drawbacks  do  not  instill  confi¬ 
dence  that  this  method  can  be  successfully  applied  to  general  fluid/structure  interaction 
problems. 
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3.3  Three-field  Approaches 

Three-field  systems,  which  utilize  two  scalar  field  variables  in  the  fluid  domain,  lie  mid¬ 
way  between  two-field  and  displacement-based  approaches  i  the  number  of  degrees  of 
freedom  present  in  the  fluid  model.  As  with  the  latter,  the  three-field  formulations  have  a 
banded  form,  which  will  usually  offset  the  disadvantages  of  having  a  larger  overall  dimension 
than  a  symmetrized  two-field  form.  The  first  three-field  system  appears  to  have  been 
proposed  by  Morand  and  Ohayon50,  who  also  noted  its  relationship  to  the  condensed  system. 
Their  approach  utilized  the  pressure  and  displacement  potential  as  fluid  field  variables,  and 
provided  a  sparse  system  in  the  form 


in  which  a  fluid  coupling  term  is  introduced,  in  this  case  of  the  form 


The  linear  dependence  of  the  three  variables  is  apparent  in  the  assembled  mass  matrix. 
Those  authors  reduced  this  system  to  a  two-field  ( S,p )  relation  in  the  form  of  equation  (32). 
Solution  difficulties  associated  with  the  singularity  of  the  assembled  stiffness  matrix  in  the 
full  system  were  not  addressed. 

Very  similar  approaches  have  more  recently  been  taken, 19,89  in  which  the  displacement 
potential  and  pressure  variables  can  be  described  independently.  This  approach  allows  the 
use  of  lower  order  interpolation  functions  for  if>,  which  is  justified  in  some  sense  since  the 
displacement  potential  is  an  integral  of  the  pressure.  As  a  consequence  of  the  possibility  of 
independent  discretization,  the  B  and  H  submatrices  need  not  be  square.  In  practice,  it 
is  probably  simpler  to  maintain  consistency  in  the  interpolation  order,  since  several  of  the 
submatrices  are  then  identical. 

In  Reference  19,  the  constant  potential  mode  solution  associated  with  the  singularity 
in  the  H  matrix,  (e)  =  c{l,l,...},  with  c  a  constant,  is  eliminated  by  forcing  the  solu¬ 
tion  eigenvectors  to  be  orthogonal  to  this  mode.  This  constraint  can  be  imposed  at  the 
variational  statement  level,  much  like  the  inclusion  of  rotational  stiffness  terms  in  the  dis¬ 
placement  variable  fluid  formulation.  The  result  is  an  added  degree  of  freedom,  A,  in  the 
assembled  system, 
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in  which 


[Be]  =  [B]{e}T  [L.\  =  [L]{e}T  (57) 

This  system  can  be  partitioned  into  two  variable  sets  (6,p,A)  and  (ip).  If  at  least  one 
constraint  in  the  ip  field  is  imposed,  then  reduced  forms  of  the  above  system,  reminiscent 
of  the  symmetrized  system  of  equation  (32)  can  be  formulated  in  only  the  first  variable  set 
for  systems  which  are  actually  singular.  This  is  an  approximation,  but  since  the  fluid  field 
variables  are  assumed  independently,  the  error  introduced  with  such  a  constraint  applied 
to  one  field  may  not  significantly  influence  the  solution  in  the  second  field.  The  reduction 
scheme  suffers  from  the  same  type  of  computational  overhead  as  the  other  symmetrization 
approaches,  since  an  explicit  inverse  is  required.  Reference  19  also  discusses  a  method 
for  solution  of  the  assembled  system  of  equation  (56),  which  would  appear  to  be  a  more 
practical  approach  in  most  cases. 

Another  three-field  system  in  which  the  (6,p,ip)  variable  set  is  utilized  was  given  in 
Reference  39.  The  assembled  system  in  this  case  is 


in  which 

(6i)<,  =  ^  VN?-VNfdn  (59) 

and  the  interpolation  functions  may  again  be  independent.  The  fluid  stiffness  submatrix 
retains  the  form 

(hi  )ij=  I  VN?  •  VNf  dCl  (60) 

Jn, 

but  may  be  of  reduced  order.  The  linear  dependence  inherent  in  this  three-field  system 
appears  in  the  assembled  stiffness  matrix,  rather  than  the  mass  matrix  as  in  equation  (54). 

While  these  three-field  approaches  provide  a  symmetric  form  for  the  coupled  problem, 
contain  fewer  fluid  variables  than  displacement  variable  formulations,  and  partially  maintain 
a  bandedness  in  the  assembled  matrices,  they  present  a  number  of  difficulties  in  terms  of 
solution  algorithm  requirements.  All  systems  have  assembled  stiffness  matrices  which  are 
singular.  In  the  case  of  equation  (58) ,  the  B  matrix  is  semidefinite,  and  this  combined  with 
the  negative  diagonal  terms  leaves  the  assembled  stiffness  matrix  negative  definite.  Also, 
as  a  result  of  utilizing  dependent  fluid  variables,  there  are  more  zero  energy  modes  in  the 
stiffness  matrix  than  can  be  removed  through  shifting. 

Despite  these  adverse  numerical  characteristics,  the  three-field  system  is  physically  well- 
posed.  Unfortunately,  current  eigensolver  algorithms  within  VAST  cannot  treat  these  sys¬ 
tems,  and  the  development  of  non-standard  techniques  will  be  required.  The  generation  of 
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the  submatrices,  based  on  isoparametric  axisymmetric  and  three-dimensiona]  elements,  and 
assembly  into  the  three-field  form  can  currently  be  carried  out  with  the  COUPLE  program. 

Static  condensation  of  equation  (58),  in  this  case  of  the  pressure  variable,  yields  the 
two-field  system 


K  +  pLQ~lLT 


-pLQ~1B1 

pBxQ-'B* 


0  111 4 1  =  1  °1 
PR i  JJ  \  0  /  1  0  / 


(61) 


which  is  essentially  in  the  form  of  the  system  of  equation  (33),  with  allowance  for  the  use 
of  mixed  interpolation.  The  comments  made  previously  on  that  system  still  apply.  As  sug¬ 
gested  in  Reference  22,  it  is  almost  always  more  efficient  to  solve  the  original  highly  banded 
system  (58)  with  a  specialized  eigensolver  than  form  and  solve  the  system  of  equation  (61). 

Another  interesting  three-field  approach  has  been  suggested  by  Olson  and  Bathe56,  in 
which  a  single  hydrostatic  pressure  variable  po  and  a  set  of  nodal  velocity  potentials  <f> 
is  used  to  describe  the  fluid.  In  this  form,  both  static  and  dynamic  problems  can  be 
solved,  and  many  fewer  fluid  variables  are  required  them  is  the  case  for  standard  three-field 
representations.  This  approach  introduces  three  new  matrices, 


(ai),  =  -^/  dfl 

P  J n. 

(62) 

r  dN+ 

M-k-si-" 
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(64) 

where  P  is  the  fluid  bulk  modulus.  The  assembled  system  becomes 
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As  was  the  result  in  the  scalar  use  of  velocity  potential,  the  effective  time  integration  of 
the  pressure  causes  the  coupling  terms  to  be  velocity  dependent.  According  to  the  authors, 
this  formulation  does  not  result  in  a  true  damping  matrix,  and  the  above  system  can  be 
solved  with  algorithms  operating  in  simple  arithmetic.  The  assembled  stiffness  matrix  is 
also  negative  definite.  The  approach  has  been  implemented  in  a  version  of  the  finite  element 
code  ADINA,48  in  which  the  determinant  search  eigenvalue  solver  system  has  been  modified 
to  obtain  solutions  of  this  problem.  The  approach  is  appealing  for  its  generality,  despite 
the  additional  matrices  required  since  these  terms  are  easily  calculated.  Unfortunately, 
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VAST  does  not  use  a  determinant  search  method,  and  further  work  would  be  required  to 
investigate  whether  modified  subspace  iteration  or  direct  iteration  methods  could  be  used 
on  this  problem.  Alternatively,  a  specialized  determinant  search  method  could  be  developed 
for  the  VAST  system. 


3.4  Modal  Coupling  Methods 


An  alternative  to  the  direct  formulations  requiring  finite  element  discretization  of  both 
structural  and  fluid  regions  is  to  treat  the  coupled  problem  using  modal  representations 
of  either  or  both  of  the  domains.  In  the  latter  case,  the  development  is  quite  analagous 
to  that  utilized  in  structural  component  modal  synthesis,  although  the  assembled  system 
for  a  coupled  fluid/structure  problem  is  usually  more  cumbersome.  The  modal  sets  re¬ 
quired  in  the  formulation  may  also  be  more  computationally  demanding  to  obtain.  Much 
of  the  modal-based  work  is  aimed  more  towards  the  accurate  prediction  of  the  acoustic  field 
response  than  the  structural  response.  Thus,  several  authors  have  used  quite  simplified 
structural  modal  representations  but  retained  a  more  detailed  fluid  model.  Regardless  of 
these  simplying  assumptions,  even  the  simplest  modal  methods  often  can  provide  a  good 
deal  of  insight  into  the  behaviour  of  the  coupled  system. 

A  simple  mixed  approach  which  retains  a  detailed  pressure-based  finite  element  fluid 
model  but  utilizes  boundary  structural  modes  was  proposed  by  Craggs.5  In  this  case,  the 
flexible  boundary  displacement  is  represented  by  a  modal  expansion  in  the  form 

u,  =  <pt$  (66) 


where  tp ,  has  an  associated  natural  frequency  o>,,  and  £  is  a  generalized  displacement.  For 
a  harmonic  forcing  function  F  at  frequency  ui,  the  displacement  response  amplitude  then 
becomes 

tf<p§F 

*  -  w2)  (67) 

and  the  boundary  velocity  amplitude  due  to  F  becomes 

•  _  .  <pJ<P$F 

*  ^Wrn(u>2  —  w2) 

These  relations  utilize  generalized  m;*ss  M  and  stiffness  K  terms  defined  by 

M  =  <pT,  M<p,  R  =  <p]  K<p,  (69) 

If  the  force  results  from  the  acoustic  backpressure  p,  then  equation  (68)  defines  a  boundary 
admittance  function  A,  where 
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This  relation  can  be  used  directly  in  the  boundary  term  of  the  standard  acoustic  functional, 
equation  (20),  to  yield  the  system 

pw*L<p,<p*  L] 

M(oj 2  —  u>2) 

This  relatively  simple  formulation  provides  a  direct  solution  to  the  coupled  problem  in 
the  limiting  cases  for  which  the  boundary  frequency  u>,  is  much  smaller  or  much  larger 
than  the  forcing  frequency,  or  equivalently,  the  fundamental  acoustic  cavity  mode.  In  the 
former,  the  boundary  is  mass  controlled,  and  the  admittance  relation  results  in  an  added 
stiffness  term,  increasing  the  acoustic  cavity  modal  frequencies.  In  the  latter,  the  response 
is  stiffness  controlled,  leading  to  an  added  mass  effect,  which  decreases  the  cavity  modal 
frequencies.  That  the  controlling  and  response  factors  seem  opposite  is  a  consequence  of 
using  the  pressure  type  of  formulation.  Solutions  for  near-resonant  cases  cannot  be  obtained 
directly  with  this  formulation. 

Resonant  boundaries  can  be  treated  with  a  more  versatile  modal  formulation,  in  which 
the  generalized  boundary  displacements  f  for  the  resonant  modes  are  incorporated  into  the 
solution  vector.  A  symmetric  system  can  be  obtained  in  the  form 
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in  which  K  is  the  generalized  resonant  modal  stiffness  defined  in  equation  (69),  and 
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represents  the  added  stiffnesses  of  all  of  the  boundary  modes  considered.  As  discussed 
in  Reference  5,  this  approach  can  treat  multiple  cavities  and  both  coupling  and  flanking 
boundaries. 

Following  a  similar  approach,  Dowell9  presented  a  genera]  form  of  coupled  modal  solution 
in  which  both  the  fluid  and  structural  boundary  motion  were  represented  with  modal  sets 
i'Pft’Pt)  respectively.  The  coupled  equations  in  this  case  utilize  non-dimensional  pressure 
expansion  coefficients  Pn  and  generalized  structural  displacements  f .  The  undamped  system 
becomes 
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in  which  the  L\  matrix  now  contains  coupling  components  defined  in  terms  of  integrated 
modal  cross-products  of  the  form 
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The  respective  modal  frequencies  in  each  domain  appear  as  diagonal  spectral  submatrices  in 
the  assembled  stiffness  matrix.  The  pressure  field  can  be  recovered  once  the  Pn  are  available. 
Following  the  methods  outlined  earlier,  this  system  can  be  symmetrized  in  several  ways. 

A  form  of  this  approach  in  which  axisymmetric  finite  element  analysis  was  used  to  define 
the  structural  and  acoustic  natural  modes  was  utilized  in  Reference  38  to  determine  the 
coupled  response  of  an  air-filled  cylinder.  As  can  often  occur  in  practice,  a  large  number 
of  structural  modes  were  predicted  in  the  vicinity  of  the  acoustic  cavity  modes,  which 
fortunately  were  widely  separated.  Although  many  of  the  structured  modes  may  have  little 
effect  on  the  coupled  response  of  the  system,  their  individual  contributions  are  unknown  a 
priori ,  and  they  must  be  maintained  in  the  problem  formulation.  Thus,  the  modal  approach 
maintains  its  effectiveness  in  the  general  case  when  only  a  very  limited  frequency  range  is 
considered. 

In  the  above  formulations,  in-vaeuo  modes  of  the  structure  have  been  utilized  with  at 
least  some  success  in  the  analysis  of  coupled  systems  which  involve  low  density  fluids.  In  fact, 
neither  the  in-vaeuo  structural  modes  nor  the  acoustic  cavity  modes  satisfy  the  true  interface 
boundary  condition  of  equation  (13),  and  this  can  lead  to  significant  error  in  problems  which 
involve  liquids.  The  work  of  Morand,  Valid  and  Oh ayon 30,34  in  particular  has  addressed 
this  issue,  in  which  mathematically  admissible  modal  basis  functions  have  been  defined 
for  the  general  problem.  For  an  arbitrary  geometry,  the  modal  sets  are  determined  via 
finite  element  analysis.  In  their  approach,  two  types  of  modal  basis  functions  are  defined: 
the  first  is  a  combination  of  the  wet  modes  of  the  structure,  where  the  fluid  is  treated  as 
incompressible,  coupled  with  quasi-static  elastic  fluid  modes;  the  second  set  is  the  acoustic 
rigid- walled  modes  of  the  cavity.  The  first  set  of  modes  satisfies  the  required  interface 
boundary  condition,  while  the  acoustic  cavity  modes  have  a  vanishing  pressure  gradient  on 
the  boundary,  and  thus  do  not  corrupt  this  interface  condition.  The  price  to  be  paid  in 
this  more  mathematically  correct  approach  is  the  increased  computation  required  to  define 
the  modified  wet  modes.  The  assembled  system  ends  up  very  much  resembling  that  in 
equation  (74) 
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in  which  the  Li,  L2  and  terms  again  represent  modal  coupling  contributions. 

A  related  approach  has  more  recently  been  taken  by  Daniels34  in  the  context  of  applying 
a  fixed-interface  component  mode  technique  to  fluid/structure  interaction  problems.  As  in 
the  work  of  Morand  and  Ohayon,  the  formulation  utilizes  modified  wet  modes  and  acoustic 
cavity  modes  to  satisfy  the  interface  boundary  condition.  This  basis  set  is  augmented  with 
up  to  two  other  types  of  modes:  a  set  of  fixed  interface  modes  (p  =  0  on  T*,)  analagous  to 
those  utilized  in  certain  formulations  of  component  mode  synthesis,  and  so-called  constraint 
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modes  in  the  fluid.  The  latter  are  incompressible  modes  which  include  the  fluid  nodal 
constraints,  and  are  included  in  the  analysis  of  fluid  domains  modified  by  the  release  of 
particular  constraints.  Again,  for  low  density  fluids,  the  in-vacuo  structural  modes  could 
be  used  in  place  of  the  wet  modes,  since  it  is  less  critical  that  the  correct  boundary  condition 
be  satisfied.  The  wet  modes  should  be  used  for  cases  involving  liquids.  The  fixed-interface 
modes  are  redundant,  but  apparently  improve  the  convergence  rate  of  the  eigensolution  of 
the  overall  system. 

With  this  seemingly  extensive  modal  set,  a  non-symmetric  modal  formulation  of  the 
coupled  problem  can  be  developed,  which  can  be  symmetrized  using  any  of  a  number  of 
the  standard  substitutions.  Because  of  its  notational  awkwardness,  the  reader  is  referred 
to  Reference  34  for  the  actual  form  of  the  assembled  system.  While  this  approach  appears 
cumbersome,  the  matrices  involved  are  generally  of  quite  small  dimension,  and  as  in  the 
other  modal  approaches,  several  of  the  submatrices  are  diagonal.  The  only  constraint  in  this 
method  appears  to  be  a  requirement  that  the  structural  stiffness  matrix  be  positive  definite. 
The  reference  contains  an  investigation  of  the  efficiency  of  the  above  form  in  the  analysis  of 
several  small  fluid/structure  problems,  comparing  the  performance  to  the  equivalent  sym¬ 
metrized  form  of  the  pressure  formulation,  equation  (34).  It  is  concluded  that  the  present 
modal  formulation  is  more  effective  from  a  purely  computational  point  of  view,  although 
as  discussed,  equation  (34)  is  no  doubt  the  most  inefficient  of  the  symmetrized  Bystems.  A 
similar  modal  approach  is  also  suggested  for  predicting  a  compressibility  correction  on  the 
pressure  field  in  a  nearly-incompressible  fluid. 

The  use  of  modal  methods  has  the  potential  advantage  of  reduced  problem  size,  which 
can  be  critical  in  many  problems  of  practical  importance.  As  well,  modal  concepts  often 
simplify  the  system  when  a  specific  frequency  range  can  be  considered.  On  the  other 
hand,  this  can  be  a  limitation  when  more  general  response  characteristics  are  required.  To 
maintain  correct  mathematical  admissibility,  the  more  demanding  calculation  of  wet  modes 
is  required,  which  will  itself  require  a  coupled  finite  element  model.  The  use  of  only  the  wet 
modes  to  represent  the  structural  response  assumes  that  the  fluid  compressibility  does  not 
affect  the  structural  mode  shape,  which  may  not  always  be  the  case.  Lastly,  as  with  the 
basic  scalar  approaches,  the  symmetrization  process  destroys  much  of  the  b finding  available 
in  the  modal  systems. 

In  light  of  these  considerations,  it  would  appear  that  the  modal  methods  will  have 
distinct  advantages  over  direct  solutions  in  problems  for  which  wet  modes  are  not  required, 
for  time  domain  solutions  which  span  limited  frequency  ranges  and  perhaps  for  modification 
type  problems.  The  direct  methods,  on  the  other  hand,  offer  an  increased  generality.  In  the 
development  of  an  effective  analysis  system  for  the  coupled  problem,  it  would  be  desirable 
to  embody  both  approaches. 
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3.5  Summary 

This  section  has  presented  four  quite  different  numerical  treatments  of  the  coupled 
fluid/structure  interaction  problem,  and  a  number  of  variations  within  each  general  form. 
To  reiterate,  two-field  scalar  formulations  in  pressure  or  displacement  potential  resulted  in 
unsymmetric  coefficient  matrices,  necessitating  secondary  symmetrizing  algorithms.  The 
symmetrization  itself  can  be  quite  inefficient,  and  each  of  the  symmetrized  systems  was 
shown  to  have  serious  limitations  in  practical  application.  Use  of  the  velocity  potential 
provided  symmetry,  but  resulted  in  the  presence  of  a  damping  matrix.  This  form  is  thus 
untractable  for  current  finite  element  procedures,  but  may  be  effective  in  time  domain  or 
direct  frequency  response  solutions.  The  scalar  analogy,  while  extending  the  usefulness  of 
conventional  finite  element  codes  to  treat  uncoupled  acoustic  cavity  problems,  suffers  from 
the  same  non-symmetric  limitation  in  its  application  to  coupled  problems.  It  is  concluded 
that  any  system  based  on  a  two-field  scalar  formulations  will  have  a  limited  applicability  for 
small  problems,  but  will  not  provide  an  effective  solution  capability  for  general  problems. 

Displacement-based  formulations  offer  bandedness  of  the  coefficient  matrices,  and  ease 
in  coupling  of  structural  and  fluid  domains,  but  suffer  from  inefficiency,  necessity  of  skew 
coordinates,  presence  of  circulation  modes  or  measures  to  eliminate  them,  and  numerical 
instability.  These  drawbacks  appear  to  preclude  the  use  of  this  approach  as  a  general 
purpose  capability. 

Three-field  systems  offer  reduced  problem  size  over  displacement-based  formulations 
and  maintain  to  some  degree  a  banded  form  in  the  coefficient  matrices,  at  the  expense 
of  adverse  numerical  characteristics  of  the  assembled  system.  These  approaches  are  easily 
formulated  from  the  basic  acoustic  and  structural  submatrices,  and  appear  to  offer  an 
effective  capability  provided  that  the  necessary  use  of  more  sophisticated  eigensolution 
algorithms  is  relatively  efficient.  It  has  not  been  possible  to  address  the  latter  point  in  this 
report. 

The  modal-based  methods  are  appealing  in  certain  circumstances,  typically  providing 
non-symmetric  systems  which  can  be  symmetrized  comparatively  easily  due  to  their  smaller 
size  and  more  diagonal  structure.  Their  lack  of  generality  could  be  a  major  drawback. 

Considering  the  alternatives,  if  a  single  approach  is  developed  for  solution  of  coupled 
problems,  one  of  the  three-field  systems  would  appear  to  be  most  appropriate,  subject  to  the 
condition  that  a  reasonably  efficient  solver  system  can  be  developed  as  well.  The  presence 
of  the  (6,#,po)  method  in  current  versions  of  ADINA  is  an  indication  that  that  approach  is 
considered  useful  and  reliable;  however,  the  eigensolution  system  required  in  that  approach 
is  quite  specialized.  If  significant  effort  need  be  directed  at  implementing  a  new  solver 
system,  it  would  be  desirable  that  it  be  as  general  as  possible.  Based  on  these  conclusions, 
in-house  and  contracted  work  is  currently  underway  on  an  original  eigensolution  capability 
for  the  (j.V'tP)  three-field  formulation.  The  complete  system  will  use  VAST  to  generate 


the  structural  submatrices,  COUPLE  for  acoustic  submatrix  generation  and  assembly  of 
the  submatrices,  and  the  eigensolver  for  solution  of  the  system.  At  the  time  of  writing, 
the  eigensolver  is  operational  for  small  systems,  but  work  remains  to  enable  large  problems 
to  be  handled  effectively.  Effort  is  being  made  to  keep  file  formats  and  details  consistent 
with  current  VAST  requirements,  and,  if  appropriate,  the  completed  eigensolver  could  be 
included  in  a  special  version  of  the  VAST  program. 

4  Example  Applications 

This  section  presents  several  example  applications  of  each  of  two  fundamental  approaches 
to  the  coupled  problem:  the  scalar  and  displacement  variable  formulations.  Unfortunately, 
only  a  limited  number  of  the  methods  discussed  in  the  previous  section  can  be  considered, 
since  many  of  them  would  require  extensive  program  development  or  modification.  Despite 
the  limited  scope,  these  simple  analyses  can  highlight  difficulties  inherent  in  the  approaches 
considered,  specifically  in  the  adaptation  or  integration  of  some  of  the  methods  for  use  with 
the  VAST  program.  Analyses  using  the  VAST-COUPLE  system  illustrate  a  symmetrized 
two-field  scalar  formulation,  and  the  VAST  program  itself  is  used  to  illustrate  the  scalar 
analogy  and  displacement-based  formulation.  The  scalar  analogy  approach  is  utilized  only 
for  prediction  of  acoustic  cavity  modes  since  it  results  in  unsymmetric  forms  of  the  assembled 
coefficient  matrices  in  the  interaction  problem. 

Although  the  COUPLE  program  can  be  used  to  construct  the  overall  mass  and  stiffness 
matrices  for  three-field  formulations,  evaluations  of  that  approach  were  not  undertaken 
at  this  time  since  testing  of  the  new  eigensolver  system  is  not  complete.  Details  of  the 
operation  of  the  COUPLE  program  are  presented  in  Appendix  A. 

4.1  Symmetrized  Pressure  Formulation 

The  VAST-COUPLE  system  has  been  used  in  the  analysis  of  several  small  coupled 
problems,  the  results  of  one  of  which  is  presented  here.  This  problem  utilizes  the  bottle 
analysed  in  Reference  44,  and  exploits  axisymmetric  modelling  both  for  the  fluid  and  struc¬ 
tural  models.  The  structural  model  consists  of  28  3-noded  axisymmetric  shell  elements, 
giving  58  nodes,  the  fluid  model  comprises  24  8-noded  axisymmetric  acoustic  elements,  giv¬ 
ing  123  nodes.  In  this  case,  all  structural  nodes  are  interface  nodes.  The  models  are  shown 
in  Figure  1.  Pressure  relief  constraints  were  applied  on  two  fluid  nodes  at  the  open  end  of 
the  bottle,  and  rigid  body  modes  of  the  bottle  were  eliminated  by  appropriate  structural 
constraints. 

Table  1  presents  the  results  for  the  first  five  vibration  modes  of  the  three  forms  of 
analyses,  with  both  air  and  water  as  the  fluid.  Modal  displacement  plots  for  four  of  the 
modes  of  the  structure  alone  and  the  coupled  structure/water  system  are  presented  in 
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Figures  2-5.  The  modes  with  air  as  the  fluid  are  indistinguishable  from  the  purely  structural 
modes,  and  are  not  presented. 
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Table  1:  Comparison  of  modal  frequency  predictions  (Hz)  for  the  uncoupled  and  coupled 
bottle/fluid  system. 

Despite  the  matching  in  the  natural  frequencies  of  the  two  individual  systems  when 
the  fluid  is  air,  there  no  apparent  interaction.  The  individual  system  modes  are  recovered 
virtually  unchanged.  This  is  not  surprising,  since  for  open-ended  systems,  the  coupling 
effect  shows  predominantly  as  an  added  mass  term,  which  is  directly  proportional  to  the 
fluid  density. 

When  the  air  is  replaced  by  water,  only  one  of  the  acoustic  cavity  modes  is  in  the 
range  of  the  lowest  structural  modes  (column  5),  but  due  to  the  higher  fluid  density,  there 
is  significant  mass-dominant  modal  coupling.  While  the  fundamental  structural  mode  is 
apparently  unaffected  by  the  presence  of  the  fluid,  it  is  suggested  that  the  particular  model 
used  here  is  not  well-suited  for  predicting  the  coupling  effects  for  this  mode.  This  is  a  result 
of  the  inconsistency  in  the  fluid  and  structural  models  at  the  closed  end.  The  introduction 
of  triangular  axisymmetric  elements  would  be  helpful  in  avoiding  this  problem.  Aside  from 
the  fundamental  mode,  the  coupled  modes  are  now  quite  distinct  from  the  purely  structural 
modes,  as  shown  in  Figures  3-5. 

Note  that  this  particular  problem  is  well  suited  for  analysis  with  the  standardization 
approach  implemented  in  the  COUPLE  program,  since  the  open  bottle  end  allows  the 
necessary  constraints  to  be  correctly  applied.  The  VAST  added  mass  model  would  provide 
a  useful  comparison  solution  for  this  problem;  however,  axisymmetric  fluid  elements  are  not 
currently  available  in  VAST. 
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4.2  Scalar  Analogy 

The  scalar  analogy  has  been  used  in  conjunction  with  the  VAST  program  for  the  predic¬ 
tion  of  the  acoustic  cavity  modes  of  a  rectangular  box.  In  this  method,  the  solid  elements 
available  in  the  VAST  program  mimic  single-degree-of-freedom  acoustic  elements.  This 
simple  analysis  highlights  some  of  the  characteristics  of  the  approach,  although  in  fact  this 
is  not  an  interaction  problem. 

The  finite  element  model  of  the  box,  which  was  also  used  as  an  evaluation  problem 
in  Reference  44,  is  shown  in  Figure  6.  In  accordance  with  the  requirements  of  the  scalar 
analogy,  the  Young’s  modulus  E  was  taken  as  1.0E+15,  Poisson’s  ratio  was  taken  as  E/2.0, 
and  a  sound  speed  in  air  of  330  m/sec  was  assumed.  The  material  density  appropriate  for 
the  analogy  method  is  thus  p  =l/e2.  In  fact,  this  density  required  scaling  upwards  to  avoid 
starting  vector  problems  in  the  VAST  program  during  the  eigenvalue  analysis,  and  certain 
actions  were  required  to  avoid  problems  with  the  unrealistic  value  of  Poisson’s  ratio.  Two 
of  the  three  degrees  of  freedom  per  node  were  constrained,  this  being  a  very  ineffective 
requirement  even  for  the  small  problem  considered  (102  constraints).  The  results  of  the 
analysis  are  shown  in  Table  2. 


Mode 

Scalar  Analogy 
(Hz) 

COUPLE+VAST 

(Hz) 

Theoretical 

(Hz) 

1 

701.8 

701.8 

699.2 

2 

1506 

1505 

1500 

3 

1541 

1541 

1398 

4 

1664 

1664 

1654 

CPU  Time  (Sec) 

232 

53 

Table  2:  Comparison  of  modal  frequency  predictions  for  the  rectangular  cavity. 

As  is  evident,  the  two  approaches  provide  comparable  accuracy,  although  the  penalty 
of  carrying  extra  degrees  of  freedom  is  reflected  in  the  increased  CPU  time  required  for 
the  scalar  analogy.  The  poor  prediction  of  the  third  mode  in  both  cases  is  a  result  of  the 
coarseness  of  the  model.44  It  must  be  noted  that  the  order  of  the  modal  predictions  of  both 
forms  of  analysis  are  sensitive  to  the  shifting  parameter  used  in  the  decomposition  process. 

4.3  Displacement-based  Formulation 

The  results  of  two  applications  of  the  displacement-based  elements  in  VAST  in  the 
solution  of  fluid  problems  are  presented  below.  The  first  is  the  prediction  of  axial  cavity 
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modes  for  the  constrained  system  shown  in  Figure  7.  To  make  the  approach  general,  the 
problem  was  modelled  with  solid  elements.  The  results,  which  are  typical  of  this  formulation, 
are  shown  in  Table  3  and  Figure  7.  True  modes  are  obtained  as  modes  2  and  8,  the  others 
being  circulation  modes,  with  frequencies  the  same  as  the  true  modal  frequencies.  A  number 
of  the  true  modes  are  not  obtained  from  the  analysis.  Similar  results  were  obtained  in 
Reference  1. 


Mode 

Frequency  (Hz) 

Characteristic 

Theoretical  (Hz) 

1 

22.45 

Circulation 

22.36 

2 

22.45 

True 

44.72 

3 

22.45 

Circulation 

67.08 

4 

22.45 

Circulation 

89.44 

5 

45.46 

Circulation 

111.80 

6 

45.46 

Circulation 

- 

7 

45.46 

Circulation 

- 

8 

39.58 

True 

- 

9 

123.30 

True 

- 

Table  3:  VAST  predictions  of  axial  natural  modes  of  the  fixed  boundary  problem  using  a 
three-dimensional  displacement-based  element. 

A  second  more  interesting  example,  taken  from  a  test  problem  used  by  Olson  and 
Bathe35  is  a  rigid  moving  boundary  problem  for  which  an  analytical  solution  can  be  ob¬ 
tained.  The  problem  schematic  and  finite  element  model  for  this  problem  are  shown  in 
Figure  8.  The  first  two  axial  modes  of  this  system  are  predicted  well  by  the  numerical 
solution,  as  shown  in  Table  4. 


Mode 

Frequency  (Hz) 

Frequency  (Hz) 

Analytical 

VAST 

1 

99.3 

99.3 

2 

1460 

1467 

In-vacuo  frequency:  1.58  Hz 


Table  4:  VAST  predictions  of  axial  natural  modes  of  the  moving  boundary  problem  using 
a  three-dimensional  displacement-based  element  (m=1.0). 
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The  analytical  values  are  solutions  of  the  transcendental  expression 


.  Ik  +  pcuiAcotf^)  ,  . 

2jr/  =  w  =  y - r-—AjU.  (77) 

where  the  relevant  parameters  are  shown  on  Figure  8.  This  closed  problem  cannot  be  solved 
using  the  symmetrized  pressure  formulation  since  the  Quid  stiffness  matrix  is  singular.  As  in 
Reference  35,  rigid  (no  bending)  vertical  motion  was  enforced  via  nodal  constraint  equations 
in  the  VAST  solution.  Although  in  this  case  skewed  coordinate  Bystems  were  not  required, 
it  is  easy  to  imagine  problems  in  which  the  boundary  conditions  would  be  very  difficult  to 
specify. 

A  modification  of  this  problem  in  which  the  axial  constraints  are  removed  at  the  end 
opposite  the  moving  boundary  can  be  used  to  test  the  capability  of  the  displacement- 
based  formulation  for  the  prediction  of  added  mass.  This  prediction  in  turn  can  be  used 
for  comparison  with  both  the  VAST  added  mass  system,33  and  the  symmetrized  pressure 
formulation  in  the  form  of  equation  32  available  in  the  COUPLE  program. 

The  spring/plate  system,  with  the  mass  m  decreased  to  0.01  in  place  of  1.0,  has  an  in 
vacuo  natural  frequency  of  15.8  Hz.  The  pressure  or  displacement  release  at  the  opposite 
end  knocks  out  the  compressibility  stiffening,  and  the  effect  will  be  of  an  added  mass.  Using 
e>  actly  the  same  model  in  each  case,  the  three  approaches  all  produce  the  same  result,  as 
shown  in  the  second  column  of  Table  5. 


Method 

Fundamental* 
Frequency  (Hz) 

Couple  +  VAST 

6.67 

VAST  added  mass 

6.68 

VAST  displacement 

6.67 

264 

iFull  end  release 
2  Partial  end  release 


Table  5:  Comparison  of  predictions  of  the  fundamental  axial  mode  for  the  partially  con¬ 
strained  moving  boundary  problem. 

While  the  agreement  in  solutions  in  the  fully  released  problem  as  described  above  is 
encouraging,  similar  results  should  be  obtained  when  a  partial  end  release  is  used.  The 
results  from  the  COUPLE  system  and  the  VAST  added  mass  systems  agree  well  and  are 
in  the  order  expected,  as  shown  in  column  3  of  the  table.  The  fundamental  displacement- 
based  solution  is  a  complicated  displacement  mode  of  high  frequency  that  does  not  model 
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the  purely  added  mass  effect  that  a  pressure  release  condition  provides.  This  points  out 
the  consideration  that  must  be  given  to  the  boundary  conditions  in  the  various  approaches; 
simple  release  conditions  may  not  be  equivalent,  as  is  indicated  here. 


5  Damped  systems 

The  neglect  of  damping  is  a  legitimate  approximation  in  many  types  of  fluid /structure 
interaction,  but  a  realistic  numerical  treatment  of  the  acoustics  of  low  density  fluids  such 
as  air  will  necessitate  a  model  capable  of  including  the  effects  of  concentrated  damping. 
Since  much  of  the  work  describe'"  in  this  document  has  been  oriented  toward  the  prediction 
of  acoustic  pressures  in  such  systems,  it  is  appropriate  to  consider  in  a  very  cursory  way 
how  damping  forces  can  be  introduced  into  finite  element  solutions  of  these  problems.  In 
all  but  a  few  special  cases,  the  inclusion  of  damping  precludes  the  use  of  the  simple  modal 
decomposition  techniques  which  have  been  the  subject  of  this  report.  In  certain  cases, 
complex  eigensolution  techniques  can  be  utilized,  but  if  damping  is  frequency  dependent, 
then  direct  frequency  response  methods  become  a  more  effective  means  to  define  system 
characteristics. 

Loss  of  acoustical  energy  in  an  acoustic  system  can  occur  in  the  form  of  reflection,  radia¬ 
tion  or  dissipation  in  absorbent  material.  Systems  which  incur  such  losses  can  sometimes  be 
treated  by  variational  methods  by  coupling  the  lossy  system  to  a  second  system  which  gains 
the  energy,  the  overall  system  remaining  conservative.  Alternatively,  a  direct  application 
of  the  weighted  residual  techniques  can  be  utilized,  these  having  particular  application  in 
cases  where  the  mean  flow  in  the  system  is  non-zero.  The  variational  approach  was  formally 
presented  by  Gladwell,31  using  what  is  commonly  called  an  adjoint  system.  The  problem 
variables  are  complex,  but  the  resulting  system  can  be  assembled  into  a  Bet  of  real  linear 
equations,  the  size  of  which  is  double  that  of  the  original  undamped  problem.  Alternatively, 
solution  algorithms  operating  in  complex  arithmetic  can  be  utilized.  A  slightly  modified 
form  of  the  adjoint  system  approach  has  been  embodied  in  the  work  of  Craggs,5-8  and 
others35,36,46,  in  the  analysis  of  reactive  and  dissipative  mufflers  and  vehicle  interiors,  and 
the  calculation  of  acoustic  impedances.  Although  the  purely  reactive  system  is  not  dissipa¬ 
tive,  the  damping  terms  have  a  similar  form  to  the  dissipative  systems,  which  are  discussed 
in  greater  detail  below. 

Two  general  methods  are  commonly  utilized  in  the  finite  element  treatment  of  absorptive 
linings.  In  the  first,  the  effect  of  the  lining  is  modeled  with  a  complex  impedance  boundary 
condition  on  the  lined  surfaces  (an  admittance  in  the  pressure-based  formulations).  This 
approach  is  attractive  in  the  sense  that  no  new  degrees  of  freedom  are  added  to  the  existing 
system,  but  is  limited  in  the  thickness  of  lining  which  can  be  modeled  accurately.  For 
thick  linings,  the  assumption  of  a  local  normal  impedance  implicit  in  this  approach  is 
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inappropriate. 


5.1  Thin  Absorptive  Linings 


In  both  variational  or  weighted  residual  formulations,  the  inclusion  of  non-homogeneous 
boundary  conditions  such  as  absorbing  boundaries  results  in  additional  boundary  integrals 
in  the  functional  to  be  minimized.  For  absorptive  linings  over  a  region  Ta>  the  integrals  are 
similar  to  the  fluid/structure  coupling  conditions  discussed  previously,  and  take  the  form 


juippv  dT 


(78) 


in  which  p  and  v  are  complex,  (i.e.  p  =  pT  -f  jpt).  As  in  the  development  presented  for 
the  modal  methods,  the  pressure  and  velocity  are  related  through  a  complex  impedance 
(admittance)  function 

*  =  AP~\p  (79) 

Substituting  this  relation  into  equation  (78)  yields  damping  terms  in  the  form 

cn  =  Apoj  f  NfNf  dT  (80) 


for  the  discretized  system.  For  purely  reactive  boundaries,  the  absorptive  linings  jure  in 
effect  replaced  by  openings  with  appropriate  characteristic  impedance  values  (typically  pc). 
An  assembled  form  of  the  damped  acoustic  system  then  takes  the  form 


[H]  -  W2[Q]  -  w[c,]  -w[Cr]  /  Pr  \  _  f  Fr  \ 

u\Cr\  [ff]  -  W2[Q]  -  «[C,j  \  p,  /  “  \  Ft  j 


(81) 


where  the  damping  matrices  are  sparse  according  to  the  configuration  of  the  lined  boundary. 
Organization  of  the  boundary  terms  into  related  groups  can  allow  the  damping  matrix  to 
be  constructed  in  a  very  compact  form.  The  standard  acoustic  stiffness  and  mass  matrices 
appear  in  equation  (81),  which  is  formulated  and  solved  for  each  value  of  frequency  to 
determine  the  system  response,  ^his  system  is  commonly  used  to  evaluate  transmission 
loss  values  for  muffler  systems. 

The  right-hand  side  of  the  assembled  system  is  now  non-zero,  and  although  in  principle 
a  standard  eigenvalue  form  of  this  system  can  be  defined26,  this  is  not  always  practical  since 
the  impedance  relation  may  be  strongly  dependent  on  frequency,  necessitating  an  iterative 
approach.  As  well,  specialized  eigensolvers  are  required,  since  the  modal  solutions  are 
complex.  The  load  terms  F  wise  from  specified  pressure  or  velocity  forcing  functions.  For 
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piston-type  problems,  where  a  boundary  nodal  velocity  vector  v  is  specified  over  a  region 
r„,  the  load  terms  take  the  form 

/,  =  pu,\C*]{Vi)  (82) 

where 

<i  =  I  *r?Nf  dr  (83) 

«'rf 

Direct  specification  of  pressure  boundary  forces  is  useful  for  modelling  hydrodynamically 
induced  loads,  such  as  occur  at  openings  of  moving  acoustic  cavities. 

An  important  requirement  in  this  approach,  and  one  which  does  not  appear  to  be 
well-discussed,  is  the  definition  of  an  appropriate  impedance  function.  Both  empirical  and 
experimental  values  have  been  utilized,  the  empirical  techniques  obviously  having  a  greater 
generality,  perhaps  at  the  expense  of  some  accuracy.  Following  Reference  12,  Craggs6,7,8 
first  defines  a  characteristic  normal  impedance  as  a  function  of  the  frequency,  gas  sound 
speed  and  the  density  and  flow  resistivity  R  of  the  damping  material.  This  characteristic 
impedance  is  then  used  to  calculate  the  impedance  of  the  lining  for  a  specific  thickness. 
The  validity  of  this  approach  is  limited  to  certain  ranges  of  frequency  (typically  200-4000 
Hz)  and  flow  resistivity  (1K-50K  Rayls).  These  ranges  are  adequate  for  most  applications 
amenable  to  finite  element  analysis. 


5.2  Thick  Absorptive  Linings 


A  second  general  approach  to  the  inclusion  of  concentrated  acoustic  damping  is  the 
direct  finite  element  discretization  of  the  absorptive  material.  This  approach  allows  the 
treatment  of  linings  of  extended  thicknesses,  which  can  behave  differently  than  thin  linings. 
In  this  case,  the  bulk  material  properties  (assumed  isotropic)  of  the  absorptive  layer  are 
imbedded  in  modified  dynamic  and  constitutive  relations  governing  the  pressure  and  velocity 
distributions  within  the  damping  material.  Within  these  elements,  these  relations  take  the 
form 

(jvpK,  +  R)u  =  -  Vp  (84) 


j<jjQ  _ 

-  —rP  —  V 

pci 


u 


(85) 


where  f I  is  the  material  porosity,  pK,  is  an  effective  material  density,  and  R  is  again  the 
flow  resistivity.  The  term  K,  is  the  structural  factor,  intended  to  account  for  the  inertia 
effect  of  the  rigid  portion  of  the  damping  material.  For  further  discussion  of  the  theoretical 
basis  of  these  concepts,  the  reader  is  referred  to  References  12  and  29. 

Applying  the  concept  of  an  adjoint  system  and  variational  formulation,  an  assembled 
system  based  on  the  relations  (84)  and  (85)  can  be  developed  for  the  dissipative  elements, 
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in  the  form 

'  [/r|-u!jr,n[<3„]  -}<*«)[<}.]  1  f  p,  \_  j  f,  \  .... 

}««J(Q.]  [/?]  -  u’K.nlQJ  j  \  p,  }  \  f,  j  ,  J 

The  Qa  matrix  is  formed  from  the  usual  functions,  with  the  effective  sound  speed  ca  used 
in  place  of  the  standard  sound  speed.  The  forcing  functions  are  defined  in  the  same  manner 
as  previously  discussed. 

For  practical  use,  these  elements  require  coupling  to  the  acoustic  elements  at  the  inter¬ 
face.  Craggs8  has  shown  that  equality  of  both  pressure  and  normal  velocity  are  required 
at  the  interface  in  order  that  these  elements  perform  correctly.  The  latter  condition  neces¬ 
sitates  a  modification  of  the  absorptive  element  matrices  of  equation  (86)  in  the  form  of  a 
multiplication  of  the  matrices  by  the  complex  factor 

j“P 

jojp  +  R 

before  assembly  into  the  complete  system.  This  factoring  and  assembly  must  be  repeated 
for  each  new  value  of  frequency. 

The  absorptive  elements  can  be  either  rigid-backed  or  allow  structural  interaction,  and 
the  other  impedance  boundaries  presented  can  be  added  to  the  complete  system.  Thus,  the 
finite  element  approach  can  be  developed  to  handle  the  major  requirements  of  a  complete 
acoustical  system  without  great  difficulty.  Certainly,  specific  features  such  as  these  damping 
models,  which  are  not  available  in  conventional  finite  element  codes  would  be  very  useful. 
Several  such  special  purpose  codes  appear  to  be  used  commercially.22*2®146 

The  development  of  a  capability  for  modelling  concentrated  damping  within  the  VAST 
program  would  be  straightforward,  although  such  a  committment  may  not  be  justified  at 
this  time.  Additional  requirements  are  essentially  limited  to  the  modification  of  existing 
algorithms  for  surface  integrations  of  interpolation  cross-products,  and  the  development  of 
effective  means  of  boundary  definition  and  overall  system  assembly. 

6  Concluding  Remarks 

This  report  has  presented  a  systematic  review  of  four  different  formulations  of  the  cou¬ 
pled  fluid/structure  eigenvalue  problem.  The  merits  and  drawbacks  of  each  general  ap¬ 
proach  have  been  discussed,  and  illustrated  with  examples  in  several  cases.  Based  on  the 
conclusions  of  this  work,  a  three-field  formulation  utilizing  structural  displacement  and  fluid 
displacement  potential  and  pressure  has  been  chosen  for  further  development.  This  work 
entails  the  generation  and  assembly  of  the  submatrices  required,  and  most  importantly  the 
development  of  an  eigensolver  capable  of  dealing  with  the  difficult  numerical  characteristics 


(87) 
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of  the  assembled  system.  At  the  time  of  writing,  a  capability  for  submatrix  generation  and 
assembly  has  been  completed,  in  conjunction  with  a  two-field  symmetrization  approach,  in 
the  program  COUPLE.  The  symmetrization  technique  is  applicable  only  to  relatively  small 
problems  and  requires  the  imposition  of  at  least  one  constraint  in  the  fluid  domain.  The 
development  of  the  eigensolver  for  the  three-field  system  is  largely  complete,  but  the  system 
remains  to  be  optimized  to  effectively  handle  larger  problems.  Future  work  will  report  on 
details  of  the  eigensolver,  and  its  application  to  some  problems  of  practical  interest.  In  the 
meantime,  the  standardization  of  the  unsymmetric  eigenvalue  problem  provides  at  least  a 
nominal  capability  for  the  analysis  of  certain  types  of  coupled  problems,  but  the  necessity 
of  matrix  inversion  limits  the  usefulness  of  this  approach. 

The  major  portion  of  the  current  work  has  concerned  itself  with  undamped  systems. 
While  the  capability  to  solve  the  real  eigenvalue  problem  is  an  important  one,  it  is  nev¬ 
ertheless  a  small  subset  of  the  general  acoustic  interaction  problem.  As  discussed,  a  more 
comprehensive  numerical  treatment  of  the  problem  will  necessarily  include  damping  models, 
and  several  methods  of  including  damping  in  finite  element  models  have  been  presented. 
Solutions  for  transient  response,  frequency  response,  complex  eigenvalue  and  nonlinear  prob¬ 
lems  are  also  necessary.  While  VAST  has  a  well-developed  capability  for  the  solution  of 
some  of  these  problems  in  the  case  of  uncoupled  systems,  much  work  remains  to  be  done  in 
the  development  of  effective  methods  of  handling  the  general  coupled  problem. 

As  a  first  step  towards  further  development,  it  is  suggested  that  capabilities  for  including 
the  effects  of  damping  be  included  in  the  solution  system.  This  is  particularly  relevant  to 
the  acoustic  fluid/structural  system  in  which  significant  concentrated  damping  can  easily 
be  introduced  in  real  problems  through  the  use  of  damping  materials.  In  addition  to  the 
concentrated  fluid  damping,  structural  damping  is  important.  This  is  often  more  of  a 
hysteritic  form  (displacement  proportional),  than  the  viscous  form  of  proportional  or  modal 
damping  currently  available. 

The  introduction  of  a  more  general  purpose  capability  to  include  either  fluid  or  struc¬ 
tural  damping  will  necessitate  development  of  algorithms  operating  in  complex  arithematic. 
While  this  is  more  costly  in  terms  of  required  computer  resources,  it  is  essential  for  the  mean¬ 
ingful  treatment  of  coupled  acoustic  fluid  /structural  interaction.  Further  effort  should  be 
directed  towards  this  goal. 
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Figure  1:  Axisymmetric  finite  element  models  of  the  bottle  and  fluid 


i 

i 

i 

MODE  NUMBER  1 


MODE  NUMBER  1 


Figure  2:  Mode  1  displacement  plots  for  the  bottle  and  coupled  bottle/water  system. 


MODE  NUMBER  2 


Figure  3:  Mode  2  displacement  plots  for  the  bottle  and  coupled  bottle/water  system. 


Figure  4:  Mode  3  displacement  plots  for  the  bottle  and  coupled  bottle/water  system. 
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Figure  5:  Mode  4  displacement  plots  for  the  bottle  and  coupled  bottle/water  system. 


Figure  6:  Finite  element  model  for  rectangular  box  of  xyz  dimension  0.236,  0.128,  0.110  m. 
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SCHEMATIC 


MODE  NUMBER  1 
2 . 245E+01  CPS 


MODE  NUMBER  2 
2 . 245E+01  CPS 


MODE  NUMBER  3 
2  245E+01  CPS 


MODE  NUMBER  7 
4.546E+C1  CPS 


MODE  NUMBER  8 
6 . 958E+01  CPS 


Figure  7:  Finite  element  schematic  and  selected  axial  modes  determined  using  a 
displacement-based  method  without  rotational  stiffness  correction. 


39 


U  U 


Figure  8:  Piston-boundary  interaction  problem  and  finite  element  model. 


40 


Figure  9:  Flow  chart  detailing  interaction  of  the  COUPLE,  VAST  and  QZ  programs  for 
eigenvalue  analysis  of  acoustic  cavity  systems. 


Appendix  A:  Program  COUPLE 


This  section  describes  the  program  COUPLE,  which  can  be  used  to  generate  the  fluid 
submatrices  required  for  a  cavity  modal  analysis  alone,  and/or  assemble  the  equivalent 
stiffness  and  mass  matrices  for  one  of  two  symmetrized  two-field  systems  or  a  three-field 
formulation  of  the  coupled  eigenproblem.  For  reference,  the  uncoupled  system  is  simply 

[\H}-^[Q}}{p}  =  {o}  (88) 

the  two  symmetrized  systems  available  are 


K  0  1  2  M  +  pLH-'L7  LH~1Q 

-  ~Q  w  -  \QH~lQ  JJ  1  p  j  -  1  0  / 


(89) 


(90) 


(91) 


Since  the  actual  application  of  the  system  (90)  has  not  been  particularity  successful  to 
date,  the  emphasis  in  the  following  discussion  will  be  on  the  formulation  of  equations  (89) 
and  (91).  The  generalized  eigenvalue  systems  of  equations  (88), (89)  and  (90)  are  solved 
using  VAST,  while  the  three-field  system  will  be  solved  with  an  eigenvalue  solution  system 
(QZ)  which  is  still  under  development.  The  characteristics  of  the  system  of  equation  (88) 
have  been  discussed  in  Reference  44,  and  will  not  be  repeated  here.  A  flow  chart  for  the 
interaction  of  the  various  programs  for  the  uncoupled  and  coupled  eigenvalue  problem  is 
shown  in  Figure  9. 

Since  structural,  fluid  and  coupled  data  files  are  all  required  at  some  stage  during  the 
formulation  of  the  coupled  system,  the  COUPLE  program  requires  three  distinct  file  prefix 
specifications.  For  the  purposes  of  this  discussion,  the  terms  PREFX,  PREFl,  and  PREF2 
denote  the  user-specified  file  prefixes  for  the  fluid  model,  structural  model  and  assembled 
model  respectively.  For  the  uncoupled  cavity  analysis,  only  the  fluid  PREFX  is  required. 

The  overall  operation  of  the  program  is  controlled  by  the  fluid  model  .USE  file,  i.e. 
PREFX. USE  in  the  current  notation.  The  format  of  this  file  is  exactly  as  required  for 
operation  of  the  VAST  program,  up  to  and  including  stiffness  matrix  modification.  The  file 
format  is  detailed  in  Reference  47.  The  key  features  of  this  file  are  the  line  of  master  control 
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switches,  which  specify  the  modules  of  the  program  to  be  executed,  and  the  additional  data 
which  control  operations  within  each  module.  Where  possible,  analogous  operations  are 
performed  by  the  COUPLE  and  VAST  program,  as  indicated  in  the  operation  summary 
given  in  Table  6. 

The  two  operations  of  fluid  submatrix  generation  and  assembly  are  quite  distinct  within 
the  COUPLE  program  and  will  be  discussed  separately.  The  structural  submatrix  genera¬ 
tion  (for  K  and  M )  is  assumed  to  have  been  performed  with  the  VAST  program,  and  will 
not  be  discussed  further  here. 


Master  contol 
switch 

COUPLE  operation 
controlled 

Module  data 
required 

1 

Element  generation 

yes 

2 

Bandwidth  minimization  (optional) 

yes 

3 

Submatrix  assembly 

yes 

4 

Constraint  imposition 

yes 

5 

Inverse  or  three-field  assembly 

no 

6 

Creation  of  .T84  and  .T85 

no 

7 

Creation  of  .T86  and  .T87 

no 

8 

Two-field  coupled  assembly 

no 

Table  6:  Master  control  switch  operations  for  the  COUPLE  program. 


Submatrix  Generation 

The  fluid  submatrix  generation  routine  is  essentially  an  extension  of  the  ACOUSTIC  pro¬ 
gram  described  in  Reference  44.  The  operations  performed  by  this  section  of  the  COUPLE 
program  are  thus: 

1.  Fluid  element  generation 

2.  Bandwidth  minimization  (optional) 

3.  Assembly  of  element  matrices 

4.  Imposition  of  fluid  constraints 

The  first  four  master  control  switches  control  execution  of  these  four  operations,  and  most 
of  the  module  control  data  in  a  standard  VAST  structural  .USE  file  is  relevant.  Restarts 
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are  possible  at  bandwidth  minimization,  assembly  and  constraint  application.  Assembly  in 
this  case  refers  to  element  assembly. 

For  fluid  submatrix  generation,  the  input  files  required  by  the  COUPLE  program  are: 

1.  PREFX.USE  -  control  file  for  submatrix  generation  and  assembly 

2.  PREFX.AMD  -  fluid  geometry  file 

The  format  of  these  files  can  be  found  in  Reference  47.  The  only  difference  occurs  in 
the  PREFX.AMD  file,  in  which  an  additional  record  containing  the  element  group  sound 
speed  must  appear  as  the  record  following  specification  of  the  element  code  and  number 
of  elements,  and  the  final  record,  which  specifies  fluid  density,  is  not  required.  The  fluid 
element  codes  are  as  specified  in  the  VAST  fluid  added  mass  section,  with  the  exception  of 
the  8-noded  axisymmetric  element.  The  element  codes  are  as  follows: 

•  IEC=3  -  8-noded  interface  element 

•  IEC=4  -  20-noded  fluid  element 

•  IEC=5  -  8-noded  axisymmetric  element 

•  IEC=8  -  3-noded  axisymmetric  interface  element 

In  all  other  aspects,  the  .AMD  file  format  is  exactly  that  required  by  the  current  VAST 
program.  For  the  uncoupled  cavity  analysis,  interface  elements  are  obviously  not  required. 

The  output  files  created  by  the  submatrix  generation  section  of  the  COUPLE  program 
depend  on  whether  the  cavity  (uncoupled)  or  coupled  analysis  is  specified  during  program 
execution.  For  cavity  modal  analysis  alone,  the  COUPLE  program  produces  only  the  assem¬ 
bled  H  and  Q  matrices.  For  the  coupled  problem,  the  coupling  L  matrix  is  also  generated, 
using  the  interface  elements  implemented  for  added  mass  computation  by  finite  element 
methods.33  The  following  files  may  thus  be  created: 

1.  PREFX.LPA  -  formatted  output  data  file 

2.  PREFX.T31  -  element  mass  matrices 

3.  PREFX.T32  -  element  stiffness  matrices 

4.  PREFX.T33  -  interface  element  matrices,  coupled  analysis 

5.  PREFX.T34  -  coupling  submatrix  L,  coupled  analysis 

6.  PREFX.T41  -  partial  geometry  file,  uncoupled  analysis 
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7.  PREFX.T45  -  bandwidth  reduction  mapping  file 

8.  PREFX.T46  -  assembled  stiffness  matrix  before  boundary  conditions  are  imposed 

9.  PREFX.T48  -  fluid  stiffness  submatrix  H 

10.  PREFX.T49  -  fluid  mass  submatrix  Q 

Note  that  the  file  extensions  are  consistent  with  those  utilized  in  the  VAST  program  for 
files  which  contain  similar  data.  All  output  files  are  unformatted  binary  with  the  exception 
of  the  .LPA  file,  which  resembles  the  .LPT  file  produced  by  VAST.  The  print  control 
parameters  in  the  geometry  and  bandwidth  control  lines  in  the  .USE  file  are  active. 

Because  constraints  are  often  required  on  the  fluid  stiffness  matrix  H ,  a  stiffness  matrix 
modification  algorithm  has  been  imbedded  in  the  COUPLE  program.  The  algorithm  is  a 
modified  form  of  that  used  in  the  VAST  program,  and  the  input,  read  from  the  ISTIFM 
section  of  the  .USE  file,  is  exactly  as  required  for  VAST.  A  minimum  of  one  constraint  is 
required  if  the  two-field  symmetrization  option  of  equation  (89)  is  to  be  utilized;  uncon¬ 
strained  systems  are  allowable  in  each  of  the  other  systems. 

Submatrix  Assembly 

The  submatrix  assembly  sections  of  the  COUPLE  program  obviously  function  quite 
differently  depending  on  the  assembly  option  specified,  but  each  option  requires  the  following 
input  data  files: 

1.  PREF1.T48  -  structural  stiffness  matrix  K  (from  VAST) 

2.  PREF1.T49  -  structural  mass  matrix  M  (from  VAST) 

3.  PREFX.T48  -  fluid  stiffness  matrix  H  (from  COUPLE) 

4.  PREFX.T49  -  fluid  mass  matrix  Q  (from  COUPLE) 

5.  PREFX.T34  -  coupling  matrix  L  (from  COUPLE) 
and  generates  the  output  files: 

1.  PREF2.T41  -  partial  geometry  control  file  for  eigenvalue  analysis 

2.  PREF2.T48  -  assembled,  coupled  stiffness  matrix 

3.  PREF2.T49  -  assembled,  coupled  mass  matrix 
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If  the  option  specified  is  assembly  of  the  three-field  system  of  equation  (91),  the  execution  is 
controlled  completely  by  the  fifth  master  control  switch  in  the  PREFX.USE  file.  Although 
VAST  is  not  being  used  for  the  eigenvalue  analysis  for  the  three-field  system  (see  Figure  9), 
these  files  maintain  the  format  currently  required  by  the  VAST  program. 

For  assembly  of  the  system  of  equation  (89),  at  least  one  constraint  must  have  been 
applied  to  the  fluid  stiffness  matrix  H .  It  is  also  recommended  that  sufficient  constraints  be 
applied  to  the  structure  to  ensure  that  the  structural  stiffness  submatrix  K  is  also  positive 
definite. 

The  computational  aspects  of  the  symmetrization  options  of  the  COUPLE  program  are 
somewhat  more  demanding  than  those  for  the  three-field  system,  and  a  number  of  inter¬ 
mediate  output  files  are  utilized,  to  both  reduce  core  memory  requirements  and  facilitate 
restarts.  With  reference  to  equations  (89)  and  (90),  the  following  files  additional  to  the 
output  files  noted  above  are  created: 

1.  PREFX.T82  -  H~x  or  M~x 

2.  PREFX.T84  -  H~XQ  or  M~XL 

3.  PREFX.T85  -  QH~lQ  or  LTM~XL 

4.  PREFX.T86  -  LH~XQ  or  KM~lL 

5.  PREFX.T87  -  pLH~xLr  or  KM~XK 

The  fifth  through  eighth  master  control  switches  in  the  PREFX.USE  file  control  the  exe¬ 
cution  of  the  symmetrization  procedure  in  these  options  within  the  COUPLE  program.  In 
this  case,  there  is  no  analogy  with  the  operations  controlled  by  these  switches  in  VAST, 
and  no  module  data  is  required  for  these  operations.  Restarts  are  possible  from  each  of  the 
control  switches. 

To  complete  the  coupled  eigenvalue  analysis  for  the  symmetrized  systems,  a  PREF2.USE 
file  will  be  required  to  control  the  decomposition  and  eigenvalue  solution  via  VAST. 

System  constraints 

The  current  version  of  the  COUPLE  program  can  assemble  systems  with  up  to  600 
structural  nodes  and  400  fluid  nodes.  This  limitation  actually  applies  only  to  the  sym¬ 
metrization  options,  since  much  of  that  system  has  been  developed  on  the  basis  of  in-core 
matrix  manipulation.  As  this  approach  is  not  well-suited  for  the  solution  of  large  problems, 
further  expansion  is  probably  unwarranted.  For  the  three-field  formulation,  the  current  size 
limitations  could  be  increased  substantially.  In  both  approaches,  it  is  important  to  exploit 
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the  bandwidth  minimization  options  in  the  program,  since  the  matrix  manipulation  algo¬ 
rithms  and  solvers  can  utilize  the  banded  nature  of  the  various  submatrices  and  assembled 
system  respectively. 

In  the  symmetrization  option  of  equation  (89),  the  inverse  matrix  is  calculated  using 
an  IMSL49  algorithm  LIN1PB  specifically  designed  for  banded  symmetric  matrices.  If 
a  diagonalized  mass  matrix  is  utilized,  a  simpler  inversion  algorithm  is  used.  Implicit 
calculation  of  the  inverse  matrices  is  not  practical  for  these  formulations. 

The  Q  matrix  is  naturally  positive  definite,  and  hence  if  the  rigid  body  modes  are 
constrained  in  the  structural  stiffness  matrix  K ,  the  assembled,  coupled  stiffness  matrix  of 
equation  (89)  will  be  positive  definite.  Thus,  frequency  shifting  will  not  be  required  in  the 
decomposition  section  of  VAST.  If  the  structure  is  free  or  partially  constrained,  shifting 
should  facilitate  solution;  however,  the  extremely  large  variation  in  the  order  of  terms  in 
the  overall  system  can  lead  to  excessive  roundoff  error,  and  the  large  band  in  the  assembled 
mass  matrix  makes  the  shifting  procedure  inefficient.  If  possible,  shifting  should  be  avoided. 
Note  that  the  presence  of  any  unconnected  fluid  nodes  will  render  the  Q  matrix  singular. 

It  should  be  noted  that  the  fluid  mass  matrix  Q  used  in  the  cavity  modal  analysis  should 
not  be  used  directly  in  the  coupled  analysis,  unless  the  reference  sound  speed  used  was  unity. 
This  in  general  is  not  possible  in  the  uncoupled  system,  since  the  entries  in  the  Q  matrix 
become  very  small,  and  the  current  version  of  VAST  will  not  generate  a  starting  vector  for 
eigenvalue  analysis.  In  the  coupled  problem,  the  assembled  mass  matrix  of  equation  (89)  is 
divided  by  p2,  which  generally  circumvents  problems  with  starting  vectors.  The  eigenvalues 
calculated  by  VAST  must  then  be  divided  by  p  to  obtain  the  true  frequencies. 

Program  input 

In  addition  to  the  options  specified  in  the  .USE  file,  a  small  amount  of  interactive  input 
is  required  by  the  program.  The  input  is  dependent  on  the  option  taken,  either  uncoupled  or 
coupled  forms  (89),  (90)  or  (91).  For  an  uncoupled  analysis,  the  following  input  is  required: 

1.  User-specified  PREFX 

2.  Uncoupled  analysis  specification 

3.  Reference  sound  speed 

For  a  coupled  analysis,  the  following  input  is  required: 

1.  User-specified  PREFX 

2.  Coupled  analysis  specification 

3.  User-specified  PREF1 
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4.  User-specified  PREF2 

5.  Fluid  density 

6.  Two  or  three-field  option  specification 

7.  H  or  M  inverse  option  (if  applicable) 

8.  Diagonal  or  consistent  mass  matrix  M  option,  (if  applicable) 

For  both  reference  sound  speed  and  fluid  density  specification,  the  program  provides  a 
choice  of  some  common  values,  or  a  user-specified  value  can  be  entered.  For  the  coupled 
analysis  specification,  the  reference  sound  speed  is  automatically  set  to  unity,  since  that 
option  cannot  be  utilized  in  the  coupled  problem.  As  previously  mentioned,  the  calculated 
eigenvalues  must  be  divided  by  the  fluid  density  to  obtain  the  true  frequencies  for  the 
coupled  analysis  when  the  form  of  equation  (89)  is  used. 

Program  Modules  and  Executable  Image  Formation 

The  current  COUPLE  program  consists  of  eleven  FORTRAN  source  code  modules  which 
contain  a  total  of  approximately  7500  lines  of  code.  The  modules  perform  operations  gener¬ 
ally  consistent  with  a  single  master  control  switch  in  the  .USE  file.  The  COUPLE  module 
controls  the  overall  program  operation.  An  executable  image  of  the  COUPLE  program  can 
be  obtained  by  compiling  and  linking  the  following  modules: 

COUPLE,  CELEMS,  ABANRD,  ASSEM,  ASTIFM,  MINVT,  MKMUL,  MLMUL, 

KMASSM,  MKASS,  KMASS3 

In  addition,  the  LINK  command  must  currently  include  the  location  of  the  object  library 
containing  the  IMSL  routine  LIN1PB,  or  a  suitable  equivalent.  The  submatrix  generation 
and  assembly  is  carried  out  in  double  precision,  but  the  matrix  manipulation  required  in 
the  symmetrization  approach  utilizes  single  precision.  The  COUPLE  program  is  currently 
compatible  with  version  5  of  the  VAST  program.  In  future,  the  symmetrization  algorithms 
may  not  be  required.  In  that  case,  the  modules 

COUPLE,  CELEMS,  ABANRD,  ASSEM,  ASTIFM,  KMASS3 
will  be  all  that  are  required.  Minor  changes  to  the  COUPLE  module  will  also  be  necessary. 
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